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1  Introduction 


Highly  flexible  aerodynamic  systems  are  ubiquitous  in  nature  and  many  engineered  applica¬ 
tions  and  there  is  increasing  evidence  that  compliance  can  confer  aerodynamic  advantages 
in  stall  and  unsteady  gusts,  and  can  extend  system  benefits  due  to  the  reduced  structural 
weight.  Furthermore,  with  increased  interest  in  highly  maneuverable  flight  vehicles,  under¬ 
standing  and  modeling  the  behavior  of  flexible  systems  during  high  angle-of-attack  and/or 
high  rate-of-change  flight  conditions  has  become  a  compelling  technical  priority. 

Aeroelastic  vibration,  or  “flutter”  has  been  studied  since  the  dawn  of  powered  flight,  and 
a  significant  body  of  research  and  understanding  has  been  assembled  [1-3].  Flutter  occurs 
when  aerodynamic  forces  induce  structural  motion.  In  many  cases,  either  aerodynamic  or 
structural  damping  is  sufficient  to  limit  these  vibrations.  However,  in  certain  situations  the 
damping  is  insufficient  and  unstable  structural  vibrations  occur,  leading  either  to  a  stable 
finite  amplitude  limit  cycle  oscillation  (LCO),  or  unstable  failure  of  the  system. 

While  the  theoretical  modeling  of  small  aerodynamic  perturbations  is  reasonably  well  un¬ 
derstood,  there  are  several  areas  that  still  demand  improved  understanding,  and  additional 
research.  These  areas  include  systems  in  which  there  is  vortex  shedding,  and  hence  signif¬ 
icant  aerodynamical  loading  from  the  free  wake.  In  these  situations,  the  modeling  of  the 
aerodynamical  loads  are  highly  nonlinear,  due  to  the  complex  evolution  of  the  vorticity  in 
the  wake.  A  second  area  where  considerable  research  is  required  is  in  which  the  aerodynamic 
structure  has  degrees  of  freedom  that  are  more  complex  than  simple  pitch  or  heave.  Exam¬ 
ples  of  this  include  membrane  wings  and  wings  operating  with  significant  motion  (such  as 
flapping  flight).  In  all  of  these  cases  there  are  few  theoretical  frameworks  to  guide  modeling 
and  design. 

One  of  the  challenges  in  developing  design  tools  and  dynamical  models  is  that  there  are 
not  many  extensive  experimental  results  with  which  to  compare.  While  there  have  been 
some  excellent  measurements  reported  on  the  behavior  of  flexible  aerodynamic  systems  (e.g. 
[4,  5])  and  several  on  unsteady  vibrations  of  flexible  MAV  wings  (e.g.  [6]),  the  drawback  of 
many  of  these  studies  that  while  they  can  serve  as  a  validation  of  computational  simulations, 
there  are  few  experiments  have  been  conducted  over  a  wide  range  of  parameters  that  isolate 
different  physical  phenomena,  and  can  be  used  to  develop  a  fundamental  understanding  of 
aerodynamic  interactions  with  flexible  structures. 

Furthermore,  analytical  and  numerical  methods  are  sorely  lacking  in  their  ability  to 
predict  vortex  growth  and  formation  due  to  large  amplitude  structural  motions.  Simple 
models  such  as  unsteady  Kutta  conditions  are  ad  hoc,  while  full  numerical  simluations  are 
prohibitively  expensive  to  run  over  a  wide  range  of  parameters.  An  intermediate,  moderate- 
order  quasi-analytical  approach  is  called  for  to  bridge  the  gap. 

1.1  Outline  of  AFOSR  program,  and  this  report 

The  two  demands  -  for  high  precision  experiments  and  good-fidelity,  moderate  cost  numer¬ 
ical  techniques  motived  our  AFOSR  program,  and  this  report  outlines  the  progress  made 
during  the  past  three  years.  In  the  next  few  pages,  we  outline  the  technical  approach  and 
achievements  in  both  experimental  and  theoretical  domains.  We  end  with  an  administrative 
summary  describing  the  various  people  who  worked  on  this  program,  and  the  conference 
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Figure  1:  Schematic  of  1-D  cyber-physical  system  (left),  along  with  the  control  system  (right). 
This  system  has  been  successfully  employed  to  explore  the  dynamies  of  systems  with  arbitrary  mass, 
stiffness  and  damping,  including  systems  with  highly  nonlinear  properties,  non-physical  properties 
and  to  experimentally  charaeterize  virtual  systems  with  the  correct  dynamics  of  those  at  a  dramat¬ 
ically  different  physical  scale  [7]. 


papers  and  presentations  that  resnlted  from  the  AFOSR  snpport.  Two  archival  jonrnal 
mannscripts  (the  hrst  in  print,  the  second  in  review)  are  attached  at  the  end  of  the  report. 

2  Experimental  approach  and  accomplishments 

2.1  Cyber-Physical  Systems 

Despite  the  interest  in  these  nonlinear  aeroelastic  systems,  extensive  stndy  is  hard  dne  to 

(i)  the  difficnlty  in  matching  the  relevant  non-dimensional  parameters  in  a  lab  environment, 

(ii)  the  expense  and  complexity  of  necessary  compntational  simulations  whose  fidelity  is 
often  uncertain,  (iii)  the  time  and  effort  necessary  to  change  the  physical  system  to  fine-tune 
its  operation  and  optimize  the  objective.  An  ideal  blend  that  combines  experimental  and 
numerical  approaches  is  the  use  of  “cyber-physical”  systems  -  the  combination  of  a  physical 
experiment  with  a  real-time  feedback  control  system  which  has  become  technically  feasible 
in  recent  years.  The  feedback  control  system  provides  a  means  by  which  a  laboratory-scale 
experiment  can  allow  for  a  detailed  experiment  that  mimics  the  dynamics  of  a  different  model 
(i.e.  one  with  different  mass,  stiffness  or  damping),  or  even  a  model  that  obeys  a  hctitious 
set  of  dynamical  constraints. 

Although  cyberphysical  systems  have  been  demonstrated  for  low-speed  system  operating 
in  water  [8,  9],  Under  the  support  of  this  program  we  designed,  constructed,  validated 
and  extensively  used  a  Cyberphysical  system  for  use  in  large  amplitude  aeroelastic  systems 
operating  in  air.  This  is  the  first  time  this  has  been  achieved.  Figure  1  shows  a  schematic  of 
the  experimental  conhguration  used  in  the  current  work,  illustrating  the  flat  plate  mounted 
to  a  virtual  spring-damper  system  in  the  wind  tunnel  test  section. 
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Figure  2:  Coefficient  of  Moment,  Cm  cls  a  function  of  angle  of  attack  (left)  and  Instability  phase 
map  (right).  The  static  coefficient  of  moment  deviates  from  the  linear  thin-airfoil  solution  leading 
to  in  intricate  series  of  bifurcations  close  to  the  traditional  static  divergence  limit. 


Rather  than  using  a  physical  torsional  spring  and  damper,  the  structural  stiffness  and 
damping  characteristics  of  the  pitching  plate  are  controlled  using  a  cyber-physical  system 
consisting  of  a  servomotor,  rotary  encoder  and  torque  sensor  coupled  to  a  real-time  digital 
control  system  (Fig.  1).  The  position,  a,  velocity,  a,  and  torque, r,  are  multiplied  by  user- 
selected  stiffness  and  damping  coefficients,  ky  and  by  respectively.  These  are  combined  with 
the  measured  torque  and  are  fed  back  to  the  servo  motor  providing  virtual  stiffness  and 
damping  and  inertia,  thus  simulating  arbitrary  structural  dynamics.  The  practical  benefits 
of  the  Cyber-Physical  system  are  chiefly  that: 

1.  The  structural  dynamics  can  be  modified  “on  the  fly”  allowing  for  rapid  and  extensive 
exploration  of  the  structural  parameter  space.  This  allows  for  fast  measurements  over 
a  wide  parameter  range. 

2.  We  can  change  structural  and  fluid  parameters  independently,  allowing  for  complete 
freedom  in  isolating  specific  physical  phenomena  and  varying  different  non-dimensional 
ratios.  For  example  the  non-dimensional  stiffness,  k*  =  2k/{pU^c^h),  can  be  kept 
constant  even  if  the  wind  speed,  Uoo,  chord,  c,  and  span,  h,  are  changed.  This  allows 
for  detailed  explorations  of  Reynolds  number  effects  while  keeping  structural  parameter 
unchanged. 

2.2  Onset  of  large-amplitude  aeroelastic  instabilities 

The  onset  of  flutter  has  been  well  documented  [1,  10]  and  in  our  work  we  have  detailed 
the  stability  boundaries  for  a  real  system  (i.e.  idealized,  but  still  with  the  full  complexities 
of  the  physical  world),  and  find  a  rich  and  detailed  map  of  the  aeroelastic  instability,  it  is 
well-known  that  the  classic  Theordorson  static  divergence  occurs  when  the  restoring  stiffness 
of  the  torsional  spring  can  no  longer  contain  the  aerodynamic  torque  as  the  plate  pitches  up. 
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However,  for  the  case  of  a  thin  flat  plate,  the  aerodynamic  coefficient  of  moment.  Cm  is  not  the 
straight  line  produced  by  Thin  Airfoil  Theory.  Rather  the  Cm{c()  curve  has  some  complexity 
(2b)  due  to  the  presence  of  a  small  separation  bubble  at  the  sharp  leading  edge,  even  at  very 
small  angles  of  incidence.  As  the  stiffness  reduces,  this  leads  to  a  rich  sub-critical  pitchfork 
bifurcation  (at  1/k*  ~  0.7)  jumping  to  the  familiar  low-amplitude  limit-cycle-oscillations 
(LCO)  that  is  familiar  from  many  classical  aeroelasticity  studies  [11-13].  For  aerodynamic 
flight  systems,  this  LCO  is  to  be  avoided.  However,  in  our  case,  we  further  reduce  the 
stiffness  (or,  alternatively,  increase  the  inverse  stiffness,  1/k*),  and  observe  a  higher-order 
Hopf  bifurcation,  leading  to  very  large  amplitude  stable  limit  cycle  oscillations,  with  pitching 
amplitudes  of  order  ±90°.  Complete  details  of  this  phase  behavior  are  documented  in  [7]. 


2.3  Pitching  kinematics  and  vorticity  dynamics 


One  key  achievement  during  the 
past  three  years  has  been  to  use 
the  Cyber-Physical  system,  in 
concert  with  high  resolution  PIV 
measurements,  to  understand  the 
large  amplitude  LCO  regime  (Re¬ 
gion  HI  in  Fig.  2).  Here,  the  plate 
oscillates  a  a  very  violent  man¬ 
ner:  at  high  frequency  (0(20)  Hz) 
and  high  amplitude  (~  ±90°). 

PIV  measurements  coupled  with 
the  torque-angle  phase  portrait 
(Fig  3)  help  to  understand  the  cy¬ 
cle. 

As  the  plate  begins  to  rise 
from  its  equilibrium  position  (£g- 
ure  3:  (1)),  the  aerodynamic 

torque  increases.  The  torque 
is  lower  than  that  predicted  by 
the  steady  moment  coefficient 
curve  [7]  due  to  the  presence  of 
dynamically-generated  vortex,  lo¬ 
cated  on  the  lower  side  of  the 
plate,  that  has  survived  from  the 
previous  cycle,  and  is  resisting  the 
positive  rotation  of  the  plate.  As 
the  angle  of  attack  increases,  the 
aerodynamic  torque  continues  to 
grow  well  beyond  the  static  stall 

angle  hnally  reaching  a  maximum  value  nearly  2.5  times  greater  than  its  static  counterpart 
(hgure  3:  (l)-(3)).  This  is  a  classic  example  of  “dynamic  stall”  wherein  the  force  and  torque 
on  the  plate  is  signihcantly  augmented  by  the  formation  and  growth  of  a  strong  LEV  on  the 


Figure  3:  Dynamics  of  large- amplitude  aeroelastic  oscilla¬ 
tions.  The  angle-torque  phase  portrait  is  shown  in  the  cen¬ 
ter  illustrating  a  characteristic  oscillation  cycle.  Each  phase 
of  the  cycle  has  a  corresponding  phase-averaged  swirl  field, 
illustrating  the  growth  and  detachment  of  the  leading-  and 
trailing-edge  vortices. 


4 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


suction  side  of  the  plate,  [14],  The  LEV  core  continues  to  grow  in  size  and  strength  as  it 
entrains  circulation  of  the  small  eddies  generated  at  the  leading  edge.  At  the  same  time,  a 
train  of  opposite-signed  vortices  is  continuously  shed  from  the  trailing  edge.  The  formation 
of  these  discrete  eddies  has  been  commonly  attributed  to  a  Kelvin-Helmholtz-like  instability 
inherent  in  the  separated  shear  layer  by  numerous  investigators  [15,  16].  The  velocity  shear 
arising  from  the  high  velocity  of  the  impinging  flow  and  the  relatively  low  velocity  of  the 
separated  free  shear  layer  at  the  sharp-edges  promotes  instability  in  the  shear  layer,  and 
therefore  leads  to  the  generation  of  the  small-scale  eddies.  Moreover,  the  shedding  frequency 
associated  with  the  small  eddies  is  much  higher  than  that  of  the  primary  vortices  which  is 
comparable  with  the  pitching  frequency. 

The  aerodynamic  torque  attains  a  maxi¬ 
mum  at  a  ~  40°  after  which  it  immediately 
begins  a  precipitous  decrease  (hgure  3- (4)). 

This  torque  decrement  is  not  due  to  vortex 
leaving  the  plate,  given  that  the  attached 
LEV  is  evident  at  a  ^  62°  (hgure  3:  (5)). 

Rather,  we  hypothesize  that  the  saturation 
of  the  LEV  circulation,  and  the  movement  of 
the  vortex  core  closer  to  the  mid-chord  point 
are  responsible  for  the  decrease  in  the  aero¬ 
dynamic  torque.  Despite  this  sharp  drop  in 
the  aerodynamic  torque,  the  plate  continues 
to  rotate  (due  to  its  inertia)  until  it  comes 
to  rest  at  a  maximum  pitch  angle,  in  this 
particular,  case,  approximately  90  degrees 
to  the  free  stream  direction.  The  shear-layer 
from  the  trailing  edge  begins  to  roll  up  and 
forms  a  trailing-edge  vortex  (TEV)  (hgure 
3:  (6)-(7)).  This  roll-up  of  the  shear  layer  is 
presumably  due  to  a  global  instability  of  the  wake  which  includes  the  inhuence  of  the  LEV 
passing  just  overhead.  Upon  reaching  a  maximum  angle  of  rotation,  we  observe  a  second 
increase  in  the  aerodynamic  torque  owing  to  the  development  of  a  secondary  vortex  in  the 
close  proximity  of  the  leading  edge  (hgure  3:  (8)-(10)).  This  secondary  vortex  is  analogous 
to  that  observed  by  [17],  who  examined  the  LEV  formation  on  a  hat  plate  undergoing  a 
translating  motion  perpendicular  to  the  incoming  how. 

2.4  Circulation  growth  and  Reynolds  number  scaling 

During  the  pitch  reversal,  a  region  of  negative  hysteresis  is  observed  due  to  the  fact  that 
the  plate’s  return  motion  is  hampered  by  the  counter  torque  induced  by  the  secondary 
vortex.  For  applications  interested  in  energy  harvesting,  the  negative  hysteresis  implies  that 
the  energy  is  being  transferred  from  the  system  into  the  wake,  and  represents  a  drop  in 
harvesting  efficiency.  We  also  note  that  the  secondary  TEV  develops  as  the  shear  layer 
emanating  from  the  trailing  edge  curls  up  in  the  wake.  Finally,  the  elastic  restoring  toque 
of  the  torsional  spring  pulls  the  plate  back  to  its  equilibrium  angle  of  attack  (hgure  3: 


Figure  4:  Universal  scaling  of  the  vortex  dynam¬ 
ics  is  achieved  using  the  local  shear  layer  velocity 
rather  than  the  free  stream.  Note  that  the  circula¬ 
tion  grows  at  a  universal  rate,  and  reaches  a  sat¬ 
urated  stage,  consistent  with  a  vortex  formation 
time  of  approximately  4- 
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(11)-(13)).  Even  at  very  small  angle,  the  flow  over  the  pressure  side  of  the  plate  does  not 
completely  reattach.  This  in  turn  hampers  the  growth  of  the  aerodynamic  torque  at  the 
beginning  of  the  second-half  of  the  cycle,  resulting  in  a  smaller  torque  value  relative  to  the 
quasi-steady  prediction.  By  symmetry,  the  pattern  repeats  during  the  second-half  of  the 
cycle. 

From  the  PIV  measurements  we  can  calculate  the  Leading  Edge  Vortex  (LEV)  growth 
and  position.  Figure  4  shows  the  growth  of  the  circulation  contained  in  the  LEV,  measured 
over  a  wide  variety  of  pitching  amplitudes,  pitch  rates  (controlled  by  varying  the  stiffness 
and  damping  of  the  structure)  and  Reynolds  numbers  (controlled  via  the  chord  and  free 
stream  velocity).  The  non-dimensionalization  is  obtained  using  the  feeding  velocity  of  the 
shear  layer  [18]: 

UsL  =  27r/c/2  4-  f/oo  sin(a(t)). 

which  takes  into  account  not  only  the  local  free  stream,  but  also  the  speed  and  direction  of 
the  leading  edge.  The  data  collapse  is  impressive,  suggesting  that  this  feeding  velocity  is 
the  appropriate  scaling  parameter.  However,  the  movement  of  the  vortex,  and  its  position 
with  respect  to  the  lifting  surface  is  critical  in  determining  the  lift,  drag  and  pitching  mo¬ 
ment  associated  with  the  LEV.  Again  the  cyberphysical  system’s  ability  to  explore  a  wide 
parameter  space  indicates  that  the  pitching  moment  is  independent  of  pitch  rate  (Fig  4b). 


2.5  Implications  for  energy  harvesting 

The  ability  to  effectively  harvest  energy  from  the  free 
stream  using  this  system  is  limited  due  to  its  single  degree 
of  freedom.  The  mechanical  energy  harvested  by  the  pitch¬ 
ing  plate  is  simply  the  area  enclosed  by  the  torque-angle 
phase  portrait  (Fig.  3),  and  this  is  easily  found  from  the 
experimental  data.  Note  that  there  is  an  area  of  positive 
energy  capture,  as  the  plate  is  pulled  up  by  the  LEV,  but 
there  is  also  a  region  of  negative  hysteresis  associated  with 
the  development  of  a  second  LEV  that  resists  the  return 
of  the  plate  to  zero  incidence  angle.  The  Cyber-Physical 
control  can  be  used  to  eliminate  this  region  of  inefficiency 
by  optimizing  the  nonlinear  stiffness  of  the  system  so  that 
the  plate  snaps  back  more  quickly  at  high  incidence  angles. 

Optimizing  the  nonlinearity  yields  a  second  torque-angle 
trace  (Fig  5)  which  eliminates  the  region  of  negative  hys¬ 
teresis  and  improves  the  power  coefficient  by  54%  over  the 
baseline  condition. 

Energy  harvesting  using  only  pitching  dynamics  is  des¬ 
tined  to  be  inefficient,  and  several  groups  have  proposed 
systems  that  include  both  pitching  and  heaving  kinemat¬ 
ics  as  efficient  mechanisms  for  energy  harvesting  (e.g.  [19- 

22]).  Key  to  the  success  of  these  techniques  is  the  ability  to  form  and  grow  an  LEV  which 
gives  additional  lifting  force  to  the  foil,  and  to  shed  it  from  the  foil  at  the  optimal  moment. 


Figure  5:  Aerodynamic  torque 
versus  angular  position  phase  por¬ 
trait  for  linear  (black)  and  cu¬ 
bic  (red)  springs.  Re=9.9- ICA . 
The  nonlinear  spring  eliminates 
the  negative  hysteresis  region  and 
increases  the  cycle  frequency,  re¬ 
sulting  in  a  54%  increase  in  the 
power  coefficent,  Cp. 
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All  of  the  schemes  proposed  by  other  groups  have  dehned  the  kinematics  and  measured  the 
energy  harvesting  efficiency.  The  example  given  here  illustrates  how  we  can  achieve  optimal 
performance  passively  using  just  the  one  degree  of  freedom,  but  the  concept  applies  equally 
well  to  systems  with  both  pitch  and  heave  motions  and  a  key  advantage  of  the  cyberphysi¬ 
cal  system  is  the  prospect  of  designing  a  purely  passive  structural  system  that  has  optimal 
energy  harvesting  characteristics  -  this  will  be  a  focus  of  the  next  phase  of  research. 

3  Theoretical  approach  and  accomplishments 

This  section  outlines  and  summarizes  the  mathematical  modeling  approach  we  pursued  to 
develop  a  framework,  which  is  extensible  beyond  the  particular  phenomenon  of  interest.  Our 
model  uses  an  approximation  of  the  Navier-Stokes  equations  inspired  by  matched  asymptotic 
methods  to  derive  a  reduced  order  model.  Because  of  the  high  Reynolds  number  of  the  flow, 
the  viscous  terms  in  the  Navier-Stokes  equations  may  be  ignored  away  from  solid  walls.  Near 
the  solid  walls,  where  the  no-slip  boundary  condition  is  imposed,  the  viscous  effects  may  not 
be  ignored  -  these  effects  dehne  a  thin  region  where  Prandtl’s  boundary  layer  approximation 
is  possible  Matching  of  the  flow  variables  between  the  two  approximations  naturally  implies 
a  vorticity  flux  shed  from  the  boundary  layer.  We  found  the  approach  effective.  However, 
the  flow  that  develops  in  the  boundary  layer  is  not  described  adequately  by  Prandtl’s  bound¬ 
ary  layer  equations;  even  moderately  high  curvatures  cause  the  approximation  to  fail.  We 
hypothesize  a  simple  £x  for  the  boundary  layer  approximation. 

This  approach  was  motivated  by  our  attempt  to  use  fluid  dynamical  models  for  the 
control  and  optimization  of  cyber-physical  systems.  Such  models  must  satisfy  the  opposing 
requirements  of  retaining  the  complexity  of  flow,  while  eliminating  the  unnecessary  details 
for  rapid  computation.  We  found  no  existing  frameworks  suitable  of  this  task. 

Numerically  solving  the  governing  Navier-Stokes  equations  using  primitive  variables  is 
prohibitively  expensive  from  the  perspective  of  flow  optimization  and  control.  The  primary 
expense  arises  from  resolving  the  boundary  layer;  its  thickness  dictates  the  number  of  grid 
points  to  be  used  for  spatial  discretization,  and  the  time  step  through  the  CFL  condition.  On 
modern  desktop  computers  one  flow  solution  in  2D  takes  at  least  hours,  and  in  3D  takes  at 
least  days.  Calculating  the  adjoint  variables  and  implementing  a  gradient  based  optimization 
algorithm  [23]  increases  the  computational  time  to  at  least  months  for  2D  problems  (years 
for  3D).  Methods  like  Large  Eddy  Simulations  model  sub-grid  scale  eddies,  but  still  have 
to  resolve  the  boundary  layer  to  accurately  predict  the  shed  vorticity,  and  hence  are  not 
exempt  from  this  estimate.  Parallel  computing  on  a  cluster  of  processors  can  reduce  the 
2D  computational  time  to  few  days  and  modern  supercomputers  may  be  able  to  reduce  the 
computational  time  for  3D  to  several  hours,  if  the  algorithms  are  effectively  parallelized,  but 
it  is  uncertain  how  data  dependencies  can  be  addressed  for  parallelization.  In  conclusion, 
optimization,  state  estimation,  and  control  are  still  out  of  range  of  computational  fluid 
dynamics. 

An  hierarchy  of  reduced  order  models  of  fluid  flow  provide  an  alternative  to  avoid  the 
computational  effort  for  the  solution  of  Navier-Stokes  equations,  at  the  expense  of  reduced 
(or  sometimes  uncertain)  accuracy.  All  these  models  identify  that  viscous  effects  in  the  flow 
may  be  neglected  except  for  the  process  of  vorticity  shedding  from  the  body.  The  simplest 
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of  these  models  [24-29]  apply  to  rigid  wings.  The  recent  advances  on  developing  the  leading 
edge  separation  parameter  [30]  also  belong  to  this  category.  They  assume  the  flow  around 
the  wings  to  be  quasi-steady  and  model  the  fluid  dynamic  force  on  them  in  terms  of  the 
instantaneous  orientation,  velocity  and  acceleration  of  the  wing.  Versions  of  such  models 
that  incorporate  unsteadiness  [31,  32]  are  also  available.  The  assumptions  underlying  these 
models  make  them  unsuitable  for  large  amplitude  unsteady  motion,  where  the  flow  can 
separate  from  multiple  points  on  the  body. 

Models  in  the  next  level  of  the  hierarchy  also  only  apply  to  bodies  with  sharp  edges, 
but  model  vorticity  transport  according  to  inviscid  dynamics  [33-40].  These  models  assume 
that  vorticity  is  shed  only  from  the  sharp  edges  of  the  body  and  then  evolves  into  a  vortex 
sheet.  The  rate  of  the  vorticity  shed  is  chosen  so  as  to  eliminate  a  singularity  of  the  flow 
at  these  sharp  edges  -  the  celebrated  Kutta  condition  from  aerodynamics.  The  shed  vortex 
sheet  is  usually  represented  by  an  array  of  discrete  point  vortices  obeying  inviscid  vortex 
dynamics  or  a  single  tightly  rolled-up  point  vortex  (governed  by  the  Brown-Michael  vortex 
dynamics).  These  methods  are  also  very  efficient.  However,  they  do  not  allow  for  the 
possibility  of  separation  at  any  point  other  than  one  of  the  sharp  edges  on  the  body.  As  a 
result,  they  cannot  be  used  to  model  vorticity  shed  from  a  smooth  body  (e.g.  an  elliptical 
wing).  Moreover,  in  practice,  vorticity  is  not  always  shed  from  all  sharp  edges  on  a  body,  for 
example,  the  flow  may  not  always  separate  from  the  leading  edge  of  a  wing.  These  methods 
cannot  a  priori  predict  the  edges  at  which  the  flow  separates.  Because  the  optimal  flows  are 
suspected  to  correspond  to  a  precise  control  of  the  instance  vorticity  shedding  switches  on 
or  off  from  one  of  the  edges,  or  of  shedding  from  another  point  on  the  surface  of  body,  these 
models  appear  to  mis-represent  or  eliminate  the  effects  underlying  the  enhanced  unsteady 
performance. 

The  most  accurate  and  computationally  intensive  methods  in  the  hierarchy  account  for 
the  vorticity  in  the  boundary  layer  as  vortex  sheets  of  variable  strength  at  or  near  the 
surface  of  the  body.  The  vortex  sheets  are  represented  as  arrays  of  point  or  blob  vortices 
and  function  to  enforce  the  no-slip  condition  [35,  41-48].  These  vortices  transport  away 
from  the  boundary  by  advection  and  diffusion,  and  leave  the  boundary  layer  where  the  flow 
separates,  thus  accurately  modeling  vorticity  shedding.  The  drawback  of  this  method  is  that 
accurate  description  of  flow  requires  computational  effort  comparable  to  direct  solution  of 
Navier-Stokes  equations.  It  is  so  because  the  distance  between  the  point  or  blob  vortices 
should  be  comparable  to  a  fraction  of  the  boundary  layer  thickness,  and  plays  a  role  analogous 
to  grid  spacing  in  direct  solution  of  Navier-Stokes  equations.  These  methods  are  not  reduced 
order  models  in  the  strict  sense,  because  they  apply  for  all  Re  and  do  not  take  advantage  of 
the  large  Re  of  the  flow. 

3.1  Our  mathematical  model 

Following  the  usual  procedure  for  matched  asymptotic  methods,  the  flow  domain  is  decom¬ 
posed  into  a  boundary  layer  of  thickness  Re”^'^^  close  to  the  solid  body,  where  viscous  effects 
are  important,  and  an  outer  region  where  they  are  negligible  (see  figure  6(a)).  The  separa¬ 
tion  of  length  scales  in  the  boundary  layer  leads  Prandtl’s  boundary  layer  equation,  which 
accounts  for  generation  and  transport  of  vorticity  in  the  boundary  layer.  As  this  vorticity 
flows  out  of  the  boundary  layer,  the  method  converts  it  to  an  array  of  point  vortices  as  a 
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Figure  6:  Schematic  showing  domain  decomposition  into  an  outer  region  and  a  boundary  layer, 
(a)  The  flow  in  the  outer  region  is  described  in  terms  of  point  vortices,  (b)  The  boundary  layer 
coordinate  systems. 


way  of  matching  the  velocity  in  the  boundary  layer  with  the  outer  flow.  The  point  vortices 
re-entering  the  boundary  layer  imparts  momentum  to  the  fluid  there,  thus  accounting  for 
re-attachment.  The  solution  automatically  accounts  for  unsteady  boundary  layer  separation, 
and  thereby  captures  the  details  of  mechanisms  such  as  dynamic  stall,  leading  edge  vortex 
shedding,  and  vortex  recapture.  The  method  has  a  fixed  computational  cost  as  a  function 
of  Reynolds  number,  presented  in  subsection  3.1.4.  This  model  was  formulated  in  2013. 


3.1.1  Outer  flow: 

The  state  of  the  flow  at  time  step  n  in  the  outer  region  is  described  using  an  array  of  point 
vortices  of  strength  Tj  at  location  Vj,  j  =  1,2,3, ... ,  N"-,  where  N'^  is  the  number  of  vortices 
at  that  time  step.  These  vortices  obey  inviscid  dynamics,  and  are  introduced  or  removed  as 
required  to  match  with  the  dynamics  in  the  boundary  layer.  The  fluid  velocity  induced  by 
these  vortices  at  this  time  step  is  given  by 


TjZ  X  {x  —  r^) 
27r\x  — 


(1) 


The  total  fluid  velocity  field  is  =  uf{x)  +  V0"'(a:),  where  the  velocity  potential  cf)^ 

satisfying  =  0,  is  added  to  satisfy  far-held  boundary  conditions  and  matching  with 

the  boundary  layer.  The  total  velocity  obeys  the  (discretized  form  of)  Euler  equations  for 
inviscid  huids 
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in  terms  of  the  vorticity  uj  =  z  •  V  x  u^,  where  the  gradient  of  pressure  imposes 
incompressibility.  From  (2),  momentum  transport  [49]  is  seen  to  arise  from  the  x  ujz 
term,  along  with  the  gradient  of  the  Bernoulli  term.  The  term  x  ujz  is  related  to  the 
hux  of  vorticity,  a  fact  used  in  the  matching  procedure. 
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3.1.2  Boundary  layer: 

The  flow  in  the  boundary  layer  is  described  using  a  body  fitted  coordinate  system.  The 
arc-length  coordinate  along  the  boundary  layer  is  denoted  6  and  the  wall  normal  coordinate 
translating  normally  with  the  body  is  rj  with  the  position  vector  x{6,ri),  as  shown  in  figure 
6(b).  The  solid  body  moves  with  velocity  Uo{6)  =  Uq{6)T' (6)  -f- Vq  (6')n”(6'),  where  and  n"' 
are  unit  basis  vectors  in  the  tangential  and  normal  directions  at  time  step  n  respectively.  The 
functions  Mq  and  Vq  describe  the  motion  of  the  body,  which  may  themselves  be  prescribed  or 
governed  by  the  dynamics  of  the  solid  body.  The  computational  method  for  calculating  flow 
around  rigidly  moving  bodies  is  described  here;  the  generalization  to  flexible  bodies  needs  to 
account  for  boundary  layer  coordinates  that  stretch;  this  detail  poses  no  technical  difficulty 
[50-52] .  The  fluid  velocity  Ubi  is  described  by  decomposing  along  the  tangential  and  normal 
direction  at  time  step  n  as  u^-y{6,vi)  =  u^{6,rj)t{6)  -f-  v^{9,ri)h{9).  The  radius  of  curvature 
of  the  boundary  is  assumed  to  be  much  larger  than  the  boundary  layer  thickness,  which 
leads  to  the  simplified  governing  equations  (the  effect  of  leading  order  body  curvature  may 
be  included  without  much  difficulty  in  this  formulation) 

Ut  +  UUg  +  VUr,  +  Pe  =  Pr,  =  0,  Ug  +  Vr,  =  0,  (3) 

where  v  =  Re“^,  and  subscripts  denote  partial  derivatives.  The  velocity  components  satisfy 
u  =  Uq  and  u  =  0  at  the  solid  wall  r/  =  0.  Asymptotic  matching  with  the  outer  region  is 
imposed  as  r/  — )■  oo.  Matching  of  normal  stress  implies  that  pressure  in  the  boundary  layer 
is  equal  to  the  pressure  in  the  outer  region,  whereas  tangential  stress  leads  to  the  shedding 
of  vorticity  as  shown  in  §3.1.3.  Vorticity  uj  =  u^  —  Vg  in  the  boundary  layer  is  approximated 
as  (u  ~  Mr;  and  its  flux  in  the  normal  direction  is  vu  —  PUr,  ~  vurj  —  A  natural  boundary 
condition  to  use  is  that  the  diffusive  flux  of  vorticity  be  zero  so  that  matching  with  the  outer 
inviscid  flow  is  possible,  i.e.  =  0  at  r;  =  Praax- 

The  boundary  layer  equation  is  discretized  in  time  as 

^n+l  _  ^  ^  ^  ^n+l  _  ^n+1  ^  ^n+1  _  ^4) 

The  viscous  term  being  the  stiffest  term  is  treated  implicitly  for  stability,  and  the  implicit 
treatment  of  the  pressure  term  is  related  to  the  outer  pressure  The  advection  terms  are 
discretized  spatially  using  upwind  differences  and  the  viscous  term  is  treated  using  central 
difference;  the  precise  form  of  this  discretization  is  chosen  so  that  matching  is  carried  out  to 
numerical  precision. 

3.1.3  Matching: 

Assuming  that  the  two  flows  agree  exactly  at  the  interface  (between  the  outer  region  and  the 
boundary  layer)  at  time  step  n,  i.e.  itQ^^(a:(6',  p^ax))  =  u'^i{9 ,  p^nux)  =  +  v^'fi,  the  exact 

numerical  matching  at  time  step  n  -|-  1  requires  careful  analysis  of  the  corresponding  terms 
in  the  equations  governing  the  flow  in  both  regions.  The  main  difficulty  is  the  consistent 
determination  of  the  unknown  pressure  to  substitute  in  (4),  because  it  depends  on  the 
outer  velocity  n  -|-  1,  which  in  turn  depends  on  the  normal  velocity  in  the  boundary  layer. 
However,  this  difficulty  is  solvable  using  a  operator  splitting  technique  presented  in  table  1. 
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Figure  7:  A  schematic  describing  the  vortex  shedding  method  in  step  9  of  table  1.  The  vorticity 
flowing  out  of  the  boundary  layer  at  point  Xk  with  flux  voj  in  time  At  through  each  element  is 
converted  to  a  point  vortex  of  strength  (top  equation).  These  point  vortices  are  placed  along 
a  curve  just  outside  the  boundary  layer  to  simulate  a  vortex  sheet  shed  from  the  boundary  layer. 
The  position  of  these  point  vortices  is  determined  by  the  distance  d^  outside  the  boundary  layer 
(middle  equation).  The  distances  dk  are  chosen  so  that  these  vortices  induce  tangential  velocity 
exactly  equal  to  the  change  in  the  boundary  layer  velocity  component  resulting  from  vorticity 
transport  to  the  outer  region  (step  7  in  table  1). 


Using  this  technique,  the  time  stepping  results  in  a  numerically  exact  agreement  between 
the  outer  and  inner  flow,  and  can  be  implemented  using  efficient  computational  techniques. 
The  outer  velocity  at  time  step  n  +  1  is  split  as 
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is  the  velocity  induced  by  the 

n+l  • 


vortices  from  time  step  n  (now  moved  to  their  new 
location),  is  the  velocity  induced  by  the  new  vortices,  and  is  the  potential  flow 


part. 

The  unknown  tangential  boundary  layer  velocity  at  time  step  n  + 1  is  also  decomposed  as 
y^n+i  _  ^n+1  ^n+1  ^n+1  ^  with  each  part  playing  a  specific  role  in  matching. 

Matching  works  by  requiring  that  each  part  of  the  outer  and  boundary  layer  (see  steps 
2-9  coded  by  backgound  color  in  table  1)  velocity  decomposition  satisfy  an  update  equation 
so  that  (i)  the  sum  of  these  update  equations  gives  (2)  and  (4)  respectively,  and  (ii)  at  the 
matching  interface  each  part  of  outer  flow  matches  with  a  part  of  the  boundary  layer  (see 
equations  in  steps  [^,  [^,  and  in  table  1  coded  by  background  color).  For  complete 
details  see  table  1. 


3.1.4  Computational  complexity: 

The  computational  effort  for  our  method  is  significantly  lower  than  direct  numerical  solution 
of  Navier-Stokes  equations,  because  only  the  boundary  layer  region  needs  to  be  discretized. 
The  number  of  grid  points  in  the  boundary  layer  region  is  independent  of  Re  because  both  the 
boundary  layer  thickness  and  grid  spacing  scale  as  The  number  of  vortices  shed  per 

time  step  is  a  small  fraction  of  grid  points  in  the  tangential  direction  because  vorticity  is  shed 
only  close  to  the  outflow  stagnation  point.  Most  of  the  computational  steps  are  broken  down 
into  computations  that  take  either  0{N)  or  O(A^logA^)  effort.  Moreover,  many  procedures 
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in  the  time-stepping  are  embarrassingly  parallelizable  so  that  extension  to  multiprocessor, 
multicore,  or  graphical  processor  unit  architectures  should  be  benehcial. 


Figure  8:  Vorticity  shed  in  an 
impulsively  started  flow  around  a 
cylinder  at  Re  =  3000.  (Left) 
The  boundary  layer  equations  (3) 
are  solved  in  a  region  of  thick¬ 
ness  bRe~^^^  around  the  cylin¬ 
der,  shown  by  the  region  be¬ 
tween  the  red  dashed  line  and 
the  blue  solid  line.  The  point 
vortices  shed  from  the  boundary 
layer  are  color  coded  according 
to  strength.  The  convective  scale 
is  used  to  non-dimensionalize 
time.  (Right)  Comparison  with 
vorticity  contour  computed  by 
Koumoutsakos  and  Leonard  [53] 
(figure  21  therein).  These  results 
were  presented  in  2015. 


3.2  Test  cases 

To  test  the  method,  we  applied  it  to  impulsively  started  flow  past  a  static  and  rotating 
cylinders,  and  stationary  ellipses  at  an  angle  of  attack  to  the  impinging  flow. 

3.2.1  Impulsively  started  cylinder  past  a  stationary  cylinder: 

A  cylinder  of  unit  radius  is  subject  to  a  flow,  which  starts  at  t  =  0  from  zero  to  a  unit 
uniform  farfleld  velocity  with  Re  =  3000.  The  resulting  flow  is  depicted  in  figure  8.  We  have 
verified  that  our  results  match  the  asymptotic  [54]  skin-friction  and  form  drag  for  small  t, 
and  the  drag  coefficient  and  vorticity  profiles  [53]  for  large  t.  The  vorticity  field  is  compared 
with  the  the  results  of  Koumoutsakos  and  Leonard  [53]  in  figure  8.  The  flow  and  the  force  on 
the  cylinder  matches  results  obtained  from  CFD  at  a  fraction  of  the  computational  expense. 
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Figure  9:  Snapshots  of  shed  vorticity  using  CFD  (left)  and  shed  point  vortices  using  our  method 
(right)  for  an  impulsively  started  uniform  flow  past  a  rotating  cylinder  for  a  dimensionless  rotation 
rate  of  2  with  Re  =  1000.  The  snapshots  agree  qualitatively,  but  quantitative  discrepancies  in  the 
patterns  formed  by  the  vorticity  are  clearly  visible. 

3.2.2  Impulsively  started  flow  past  ellipses  at  angle  of  attack  and  rotating  cylinders: 

While  the  impulsively  started  flow  past  a  stationary  cylinder  agreed  well  with  direct  nu¬ 
merical  solution  of  the  Navier-Stokes  equations,  the  agreement  is  due  to  a  special  condition 
satisfied  in  this  case  but  not  shared  by  many  other  common  cases  of  interest.  The  impul¬ 
sively  started  flow  past  ellipses  oriented  at  an  angle  of  attack  and  rotating  cylinders  represent 
cases  which  show  disagreement  between  results  of  our  model  and  flow  computed  using  direct 
solution  of  Navier-Stokes.  We  systematically  investigate  the  source  of  the  disagreement  and 
determine  the  modification  to  Prandtl’s  boundary  layer  required  for  minimizing  it. 

First  consider  an  impulsively  started  flow  past  a  rotating  cylinder.  The  flow  computed 
using  CFD  agrees  well  with  the  flow  computed  using  our  model  for  slow  rotation  rates.  How¬ 
ever,  for  faster  rotation  rates,  the  results  from  our  model  disagree  with  the  direct  numerical 
solution  of  Navier-Stokes.  Figure  9  shows  snapshots  of  vorticity  computed  using  CFD  and 
distribution  of  point  vortices  computed  using  our  model  for  Re  =  1000  and  a  rotation  rate 
of  2,  which  shows  a  case  on  the  verge  of  this  disagreement. 

Since  the  only  possible  inaccurate  component  of  the  model  is  the  boundary  layer  approx¬ 
imation  (the  vortex  dynamics  accurately  represent  Navier-Stokes  equations  in  the  region 
where  viscosity  is  negligible),  the  only  possible  cause  for  the  disagreement  is  the  boundary 
layer  approximation.  We  have  separately  verified  from  the  early  evolution  of  the  flow,  when 
the  outer  region  is  irrotational  and  therefore  devoid  of  vorticity,  that  the  evaluation  of  poten¬ 
tial  flow  in  the  outer  region  forced  by  the  boundary  layer  dynamics  is  accurate.  Therefore, 
the  only  remaining  possibility  of  the  boundary  layer  approximation  being  inaccurate  must 
be  true. 

A  careful  inspection  of  the  boundary  layer  approximation,  coupled  with  numerical  ex¬ 
periments  examining  the  failure  of  our  model  allows  us  to  improve  upon  the  boundary  layer 
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Figure  10:  Snapshots  of  vortic- 
ity  field  computed  using  CFD  (left) 
and  point  vortex  distribution  com¬ 
puted  using  our  model  for  impul¬ 
sively  started  flow  past  an  ellipse. 


approximation.  Within  the  boundary  layer  approximation,  the  separation  of  length  scale 
holds,  and  so  does  the  neglect  of  the  axial  viscous  stress  term  and  the  dominance  of  the 
tangential  over  normal  velocity  component.  The  least  accurate  approximation  is  the  ne¬ 
glect  of  the  pressure  variation  along  the  thickness  of  the  boundary  layer;  the  contribution  of 
the  centripetal  acceleration  to  the  normal  pressure  gradient  may  lead  to  noticeable  errors. 
Thus  we  propose  to  modify  the  normal  momentum  balance  equation  the  boundary  layer 
approximation  as 


Ut  +  UUg  +  VUrj  +  Pe  =  l^Urirj,  Prj  =  Ug  +  =  0.  (5) 

As  an  example  consider  the  flow  around  an  ellipse  of  aspect  ratio  e  at  an  angle  of  attack 
to  the  impinging  flow.  The  maximum  curvature  on  the  ellipse  surface  scales  as  1/e  along 
the  major  axis,  where  the  maximum  velocity  of  0(l/e^/^).  The  magnitude  of  the  normal 
pressure  gradient  scales  as  1/e^,  and  therefore  the  difference  in  pressure  across  the  boundary 
layer  thickness  scales  as  0(e“^Re“^'^^).  Thus  an  aspect  ratio  as  much  as  Re^'^^  can  violate 
the  boundary  layer  approximation. 

To  test  this  hypothesis,  we  implemented  our  model  for  the  case  of  an  impulsively  started 
ellipse.  The  point  vortex  distribution  computed  using  our  model  with  the  boundary  layer 
approximated  by  (3),  is  shown  in  hgure  10,  and  is  compared  with  results  obtained  from  CFD. 
We  have  not  yet  implemented  the  modified  boundary  layer  model  (5),  but  comparison  of 
flows  computed  using  (3)  and  (5)  provides  a  systematic  way  of  comparing  the  effectiveness 
of  our  hypothesis.  With  this  modification,  the  method  retains  the  sparse  structure  of  the 
model  that  facilitates  rapid  computation. 
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4  Summary  of  accomplishments 

Understanding  and  optimizing  coupled  fluid-structure  systems  for  different  applications  is 
challenging.  Experiments  capture  the  real  physics,  but  are  not  easily  employed  over  a  wide 
range  of  scales,  speeds,  stiffness  and  damping.  Computational  efforts  can  “dial  in”  arbitrary 
parameter  ranges,  but  are  time-consuming  and  often  don’t  capture  realistic  physics,  often 
due  to  Re-effects.  The  key  accomplishments  achieved  during  the  past  three  years  of  AFOSR 
funding  are: 

Implementation  of  Cyber-Physical  system.  We  have  implemented  the  Cyber-Physical 
system  in  air  (previous  implementations  have  been  in  water)  operating  at  high  fre¬ 
quency  and  with  high  accuracy.  The  dynamics  of  the  system  are  well-characterized 
and  robust. 

Identification  of  Instability  boundaries.  The  phase  map  delineating  stable  behavior  as 
well  as  low-  and  high-amplitude  limit  cycle  oscillations  (LCOs)  has  been  mapped  out 
and  carefully  characterized. 

Physics  of  high  amplitude  oscillations.  The  high-amplitude  LCO  regime  has  been  care¬ 
fully  measured  and  each  phase  associated  with  the  formation,  growth  and  detachment 
of  the  leading  edge  vortex  associated  with  the  separation  at  the  sharp  leading  edge. 

Leading  edge  vortex  growth  and  scaling.  The  dynamics  of  vortex  growth  have  been 
characterized  and  scaled  over  a  wide  range  of  pitching  behavior  and  a  moderate  range 
of  Reynolds  numbers.  Although  this  is  not  complete,  a  consistent  description  has  been 
identified,  but  still  needs  to  be  rigorously  validated. 

Development  of  new  model  for  high-Re  flows.  Our  model  uses  an  approximation  of 
the  Navier-Stokes  equations  inspired  by  matched  asymptotic  methods  to  derive  a  re¬ 
duced  order  model.  Because  of  the  high  Reynolds  number  of  the  flow,  the  flow  is 
assumed  inviscid  away  from  the  walls,  but  near  the  walls,  where  the  no-slip  boundary 
condition  is  imposed,  the  viscous  effects  give  rise  to  a  thin  boundary  layer.  Matching 
of  the  flow  variables  between  the  two  approximations  naturally  implies  a  vorticity  flux 
shed  from  the  boundary  layer. 

4.1  Remaining  questions  and  future  work 

Despite  these  successes,  several  questions  remain,  and  the  Cyber-Physical  approach  provides 
a  unique  opportunity  to  explore  a  wide  range  of  unsteady  aerodynamics  and  vortex-structure 
interactions.  The  key  research  questions  that  we  identify  are: 

Scaling,  the  Cyber-Physical  system  allows  us  to  separate  structural  and  fluid  parameters 
with  unprecedented  flexibility  (for  an  experimental  system).  Some  indications  on  scal¬ 
ing  have  been  identified,  but  more  needs  to  be  done.  What  are  the  relevant  roles  of 
pitching  amplitude,  rate  and  acceleration  as  well  as  Reynolds  number,  and  structural 
inertia? 
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Three-dimensional  effects  Current  experiments  have  been  deliberately  restricted  to  two- 
dimensional  phenomena  that  emphasize  the  LEV.  What  is  the  role  of  three-dimensional 
effects,  particularly  spanwise  flows  that  are  present  in  swept,  rotating  and  flapping 
wings? 

Energy  harvesting.  How  can  a  passive  structural  characteristics  be  used  to  replace  pro¬ 
scribed  kinematic  motion,  and  how  can  we  optimize  the  vortex  dynamics  for  energy 
capture. 

Model  development.  The  flow  that  develops  in  the  boundary  layer  is  not  described  ade¬ 
quately  by  Prandtl’s  boundary  layer  equations;  even  moderately  high  curvatures  cause 
the  approximation  to  fail.  We  hypothesize  a  simple  modihcation  for  the  boundary  layer 
approximation  by  taking  into  account  the  component  of  centripetal  acceleration  in  the 
normal  direction,  which  can  be  implemented  at  minimum  additional  computational 
expense. 


16 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


,  Advect  vortices  (fc  =  1,  2, . . . ,  A"")  according  to  4.  Compute  satisfying  (a;' 
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Table  1:  Operator  splitting  time- stepping  algorithm.  The  veloeity  and  pressure  are  split  in  parts  as  shown  on  the  first  line,  and  the 
steps  1-12  are  sequentially  earried  out.  Left  eolumn  denotes  operations  on  the  outer  region,  and  right  on  boundary  layer  variables.  The 
equations  governing  eaeh  of  these  parts  are  shown  in  Azure  and  Almond  color  backgrounds,  and  the  relations  between  the  boundary 
layer  and  outer  variables  are  highlighted  in  other  color  backgrounds.  The  last  line  summarizes  the  matching  of  the  velocity  components 
at  the  interface  g  =  rjrnax- 
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The  dynamics  of  an  elastically  mounted  flat  plate  in  a  uniform  stream  undergoing  flow- 
induced  pitching  oscillations  is  studied  with  the  aid  of  a  cyber-physical  system  which  is 
capable  of  simulating  arbitrary  structural  dynamics.  The  system  measures  the  angular 
position  and  velocity,  and  combines  them  with  appropriate  gains  to  provide  a  torque  to  a 
computer-controlled  servomotor  that  emulates  arbitrary  torsional  stiffness  and  damping. 
A  series  of  experiments  were  carried  out  over  a  wide  range  of  parameter  space  and  the 
results  demonstrate  that  inertial  scaling  of  the  stiffness  and  damping  effectively  captures 
the  system  behavior.  As  the  torsional  stiffness  decreases,  the  system  exhibits  several 
bifurcations  in  its  dynamical  behavior.  Firstly,  a  subcritical  transition  through  a  saddle- 
node  bifurcation  results  in  small-amplitude  asymmetric  limit-cycle  oscillations  and  later, 
a  subcritical  transition  through  a  Hopf  bifurcation  gives  rise  to  large-amplitude  symmetric 
limit-cycle  oscillations.  At  very  low  stiffness  the  plate  does  not  oscillate.  Energy  harvesting 
from  the  flow  is  only  ~  1%  efficient,  but  is  found  to  be  optimized  at  a  non-dimensional 
stiffness  close  to  unity  while  increasing  the  Reynolds  number  is  found  to  extend  the  range 
of  damping  over  which  large  scale  oscillations  and  energy  harvesting  can  be  sustained. 
Although  velocity  measurements  are  not  made,  the  details  of  torque-position  phase  plane 
are  used  to  infer  the  details  of  the  formation  time  and  stability  of  the  leading-edge  vortex 
associated  with  the  rapidly  pitching  plate. 

©  2015  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

The  unsteady  motion  of  an  elastically  supported  structure  immersed  in  a  fluid  flow  has  been  a  problem  of  technical  and 
scientific  interest  for  decades,  originating  from  the  aeroelastic  excitation  of  aircraft  wings,  marine  structures,  bridges,  and 
even  traffic  signs  (Bisplinghoff  et  al.,  1957;  Blevins,  1990).  In  recent  years,  harvesting  of  the  aeroelastic  motions  has  also  been 
explored  for  energy  conversion  applications  (Bernitsas  et  al.,  2008;  Frayne,  2011 ),  The  onset  of  small  amplitude  motions  can 
be  predicted  using  classical  theories  (e.g.,  Theodorsen,  1935)  for  which  linear  potential  flow  assumptions  remain  valid. 
However,  for  larger-amplitude  motion,  vortex  growth  and  separation  dominate  the  flow  both  on  the  body  and  in  the  wake, 
and  the  linear  models  are  unable  to  predict  the  details  either  of  the  vortex  dynamics  or  the  unsteady  forces  that  the  flow 
imparts  to  the  elastic  structure.  Developing  a  more  sophisticated  model  is  challenging,  and  part  of  the  difficulty  lies  in 
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Nomenclature 

f* 

dimensionless  mass  moment  of  inertia 

b* 

dimensionless  linear  damping 

H„o 

free-stream  velocity  (m/s) 

k* 

dimensionless  torsional  spring  stiffness 

c 

chord  (m) 

Re 

Reynolds  number 

h 

span (m) 

damping  ratio 

s 

plate  thickness  (m) 

OJjj 

undamped  natural  frequency  (rad/s) 

P 

kinematic  viscosity  of  the  fluid  (m^/s) 

(Op 

pitching  frequency  (rad/s) 

p 

fluid  density  (kg/m^) 

reduced  pitching  frequency 

t 

time  (s) 

damped  natural  frequency  (rad/s) 

a 

angle  of  attack  (rad) 

reduced  damped  natural  frequency 

«o 

equilibrium  angle  of  attack  (rad) 

T 

period  of  oscillation  (s)  T  =  Inlcop 

^max 

maximum  angular  position  (rad) 

aerodynamic  torque  (N  m) 

a 

time-averaged  mean  angle  of  attack 

W 

average  work  per  cycle  Q) 

Aa 

peak  amplitude  (rad)  Aa=  Umax -a 

N 

number  of  cycles 

I 

mass  moment  of  inertia  (kg  m^) 

Cp 

power  coefficient 

k 

torsional  spring  stiffness  (N  m/rad) 

Cm 

steady  moment  coefficient  about  the  mid- 

b 

linear  damping  coefficient  (N  m  s/rad) 

chord  point 

finding  a  unifying  scaling  of  the  leading-edge  vortex  (LEV)  formation  time,  strength  and  stability  that  can  be  generalized  to 
various  geometries  and  kinematics.  Our  current  ability  to  model  large-amplitude  oscillations  due  to  fluid-structure 
instabilities  remains  insufficient,  albeit  a  number  of  experimental  studies  reported  in  the  archival  literature  (e.g.,  Dimitriadis 
and  Li,  2009;  Razak  et  al.,  2011;  Amandolese  et  al.,  2013).  Significant  efforts  have  been  made  to  elucidate  the  mechanism  by 
which  a  classical  linear  flutter  transforms  into  a  complex  bifurcation  behavior  that  leads  to  violent  oscillations  induced  by 
the  dynamic  stall  phenomenon.  Although  the  physical  mechanism  of  the  stall  flutter  is  well  understood,  it  still  poses  a  great 
challenge  to  accurately  correlate  the  aerodynamic  loads  due  to  unsteady  vortex  formation  and  shedding  with  the  resulting 
dynamics  of  an  elastically  supported  structure. 

Recent  studies  that  focused  exclusively  on  characterization  of  the  LEV  formation  have  used  predefined  kinematic  motion  of  a 
sharp-edged  flat  plate  to  explore  the  forces  induced  by  pitching  (Granlund  et  al.,  2013;  Jantzen  et  al.,  2014),  rotation  (DeVoria  and 
Ringuette,  2012;  Wojcik  and  Buchholz,  2014),  and  plunging  (Rival  et  al,  2009;  Ford  and  Babinslcy,  2013).  Numerous  studies  on 
pitching  and/or  plunging  wings  with  sufficient  forcing  have  identified  reduced  frequency  as  the  primary  scaling  parameter  that 
determines  the  flow  evolution  and  unsteady  force  (see,  e.g.,  Baik  et  al,  2012;  Jantzen  et  al,  2014).  However,  the  measurements  of 
Granlund  et  al.  (2013)  have  shown  that  the  pitch  rate  contribution  to  unsteady  force  is  largely  dependent  on  the  location  of  the 
pitch  axis,  presumably  because  the  vorticity  distribution  and  non-circulatory  loads  on  the  wing  vary  sensitively  with  the  pitch 
axis,  as  has  been  demonstrated  by  Yu  and  Bernal  (2013).  Therefore,  the  applicability  of  scaling  relations  based  solely  on  reduced 
frequency  (such  as  those  of  01  et  al,  2010;  jantzen  et  al,  2014)  may  be  quite  limited.  More  detailed  characterization  of  large- 
amplitude  pitching  plate  is  still  required  to  determine  a  more  robust  scaling  parameter  that  captures  the  essential  features  of  the 
unsteady  aerodynamics  with  variations  in  plate  configurations  and  kinematics. 

It  is  nevertheless  well-established  by  now  that  the  LEV  formation  on  an  unsteady  flat  plate  is  responsible  for  significant 
enhancement  in  lifting  force  that  far  exceeds  its  static  counterpart.  Hence  from  the  energy  harvesting  point  of  view,  one  may 
surmise  that  the  LEV  dynamics  can  also  have  a  significant  impact  on  the  amount  of  fluid  energy  that  can  be  harvested  by  the 
elastic  structures  through  large-scale  coupled-mode  flutter  mechanisms.  The  numerical  investigations  of  Peng  and  Zhu  (2009) 
and  Zhu  and  Peng  (2009)  reinforced  the  importance  of  vorticity  control  mechanism,  wherein  the  energy  of  LEVs  is  partially 
recovered  through  favorable  vortex-body  interactions,  in  enhancing  the  system's  overall  energy  harvesting  capacity. 

Experimental  studies  of  energy  harvesting  are  less  common  in  the  archival  literature.  Part  of  the  difficulty  lies  in  creating 
a  physical  experiment  in  the  laboratory  scale  that  possesses  the  characteristics  for  optimal  energy  harvesting.  One  approach 
to  overcome  this  issue  is  the  use  of  a  “cyber-physical”  system  in  which  the  mechanical  properties  of  the  system 
(e.g.  stiffness,  damping)  are  generated  using  a  control  algorithm  that  measures  the  state  of  the  system  and  feeds  back 
restoring  torques  and  forces  so  as  to  simulate  the  appropriate  structural  stiffness  and  damping.  Thus  far,  the  cyber-physical 
systems  reported  have  been  used  to  investigate  vortex-induced  vibration  (VIV)  behavior  of  elastically  mounted  marine 
cables  and  cylinders  in  water  (Hover  et  al,  1997;  Mackowski  and  Williamson,  2011,  2013;  Lee  et  al,  2011 ;  Lee  and  Bernitsas, 
2011).  The  implementation  in  water  typically  allows  for  testing  in  relatively  low  speed  flows  (and  hence  requiring  low 
bandwidth  control  systems)  while  still  being  able  to  explore  moderate  Reynolds  number  regimes  that  range  from  0(10^)  to 
0(10^).  With  the  exception  of  the  preliminary  results  presented  by  Song  (2013),  cyber-physical  implementations  of  pitching 
and/or  plunging  airfoils  or  of  systems  using  air  as  the  working  fluid  have  not  yet  been  reported. 

In  the  current  paper,  we  take  advantage  of  the  wide  range  of  parameters  that  can  be  explored  using  a  cyber-physical 
system  to  report  on  the  flow-induced  oscillations  of  a  flat  plate,  mounted  at  mid-chord  and  immersed  in  an  airstream.  This 
geometry  has  the  advantage  of  simplicity;  the  sharp  leading  and  trailing  edges  fix  the  separation  points  while  the  single 
degree  of  freedom  means  that  the  available  parameter  space  (stiffness,  damping  and  Reynolds  number)  is  constrained  so 
that  we  can  fully  understand  the  fluid-structure  interaction.  The  objectives  of  the  present  investigation  are  three-fold;  (1 )  to 
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axis  of  rotation  \  \ 


Fig.  1.  A  schematic  of  the  oscillating  flat  plate  (with  1  DOF).  The  axis  of  rotation  is  placed  at  the  mid-chord  point  to  ensure  symmetry,  and  the  sharp  edges 
of  the  plate  fix  the  flow  separation  points  at  the  leading  and  trailing  edges. 

demonstrate  the  successful  implementation  of  a  high-bandwidth  cyber-physical  system  in  a  wind  tunnel,  and  confirm  its 
capability  to  simulate  a  wide  range  of  structural  stiffness  and  damping,  (2)  to  determine  the  scaling  of  the  aeroelastic 
behavior  with  respect  to  the  geometric,  elastic  and  fluid  parameters,  and  (3)  to  assess  the  system's  ability  to  harvest  energy 
from  a  uniform  airflow  as  a  function  of  these  parameters. 


2.  Experimental  system  and  dimensional  scaling 

2.1.  Schematic  description  and  dimension  scaling 


Fig.  1  illustrates  a  flat  plate,  with  chord  c,  span,  h,  thickness,  S,  and  mass  moment  of  inertia,  I,  free  to  rotate  in  pitch  about 
its  mid-chord  in  a  uniform  airflow  of  density  p  and  magnitude  U^o-  The  structural  behavior  of  the  plate  is  determined  by  a 
torsional  spring  and  damper,  characterized  by  coefficients  k  and  b  respectively.  The  dynamics  of  this  single-degree-of- 
freedom  plate  in  response  to  aerodynamic  forcing,  rp  can  be  modeled  as  a  forced  torsional  harmonic  oscillator 

la(t)  +  ba(t)  +  k{a(t)-ao)  =  Tf{t)  (1) 


where  a  is  the  rotational  position  of  the  plate,  Og  is  the  equilibrium  angle  of  the  torsional  spring  and  dots  represent  time  derivatives. 

In  this  study,  we  are  primarily  interested  in  very  strong  fluid-structure  interactions  which  drive  large-scale  structural 
motions.  For  such  interactions,  we  hypothesize  that  the  dominant  scaling  force  will  be  related  to  the  inertial  forcing  by  the 
fluid,  (Shiels,  1998;  Gharib,  1999).  With  this  assumption,  the  reduced  torsional  spring  stiffness,  k*,  which  gives  the  ratio 
between  the  elastic  restoring  force  and  aerodynamic  force  can  be  written  as 


k*  = 


2k 

pUic^h 


(2) 


Similarly,  the  dimensionless  damping,  b*,  which  describes  the  ratio  between  the  structural  damping  force  and  aerodynamic 
force  is  given  by 


b*  = 


2b 

pU^c^h’ 


(3) 


and  the  dimensionless  inertia,  f  (often  called  the  “mass  parameter”)  which  compares  the  plate  inertia  to  that  of  the 
surrounding  fluid  is  given  by 


I* 


21 

ptrili 


(4) 


Alternative  scalings  have  been  proposed,  such  as  the  traditional  reduced  velocity  and  damping  ratio,  that  use  purely  structural 
parameters  (i.e.  the  structural  mass  and  the  natural  frequency  time  scale)  as  utilized,  for  example,  by  Khalak  and  Williamson 
(1999).  However,  such  formulation  lacks  a  clear  connection  to  the  fluid  dynamical  aspects  of  the  problem  and  causes  the 
normalized  structural  parameters  to  be  coupled.  Lastly,  provided  that  the  Reynolds  number  is  large  (in  the  present  study,  it  is 
0(10^)),  and  the  separation  points  are  well-defined  (ensured  in  the  current  study  by  the  use  of  a  sharp-edged  flat  plate),  it  is 
reasonable  to  assume  that  the  viscosity  plays  only  a  minor  role  in  the  overall  dynamics.  For  these  reasons  a  viscous  time  scale  is 
probably  not  an  appropriate  choice.  The  validity  of  the  current  normalization  will  be  tested  in  the  accompanying  experiments. 
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2.2.  Experimental  implementation 

Fig.  2  shows  a  schematic  of  the  experimental  configuration  used  in  the  current  work,  illustrating  the  fiat  plate  mounted 
to  a  virtual  spring-damper  system  in  the  wind  tunnel  test  section.  The  closed-return  wind  tunnel  has  a  test  section  of 
0.61  m  in  width,  0.61  m  in  depth,  and  1.22  m  in  length.  The  turbulence  level  was  previously  measured  to  be  less  than  0.1%  of 
the  free-stream  velocity  (Song  et  al.,  2008).  Titanium  plates,  1  mm  In  thickness,  were  attached  to  a  steel  shaft  with  a 
diameter  of  1  cm  mounted  vertically  between  the  floor  and  ceiling  of  the  test  section.  Two  different  configurations  were 
tested,  with  chord  dimensions  of  0.1  m  and  0.15  m.  In  both  cases,  the  plate  span  was  0.53  m,  giving  aspect  ratios  of  5.3  and 
3.5.  The  experiments  were  performed  at  free-stream  velocities  ranging  from  10  m/s  to  25  m/s,  corresponding  to  Reynolds 
numbers,  based  on  the  chord  length,  ranging  from  0(10'*)  to  0(10^).  Since  the  chord  represented  no  more  than  25%  of  the 
test  section  width,  blockage  and  side  wall  boundary  layer  effects  were  deemed  insignificant. 
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Fig.  2.  Experimental  setup:  a  computer-controlled  servomotor  rotates  the  shaft-plate  assembly  by  applying  a  specific  amount  of  torque  determined 
numerically  to  realize  desired  structural  dynamics. 


2.3.  Cyber-physical  system  control  scheme 

Rather  than  using  a  physical  torsional  spring  and  damper,  the  structural  stiffness  and  damping  characteristics  of  the 
pitching  plate  were  controlled  using  a  cyber-physical  system  consisting  of  a  servomotor,  rotary  encoder  and  torque  sensor 
coupled  to  a  real-time  digital  control  system  (Fig.  3).  The  top  end  of  the  shaft  was  coupled  to  a  brushless  servomotor 
(SM233AE,  Parker  Compumotor,  Rohnert  Park,  CA)  via  a  reaction  torque  sensor  (TFF400  Futek,  Irvine,  CA)  while  the  bottom 
end  of  the  shaft  was  attached  to  a  rotary  encoder  (HB6M  Optical  Encoder,  US  Digital,  Vancouver,  WA).  A  tachometer  (ETCH2, 
USDlgltal,  Vancouver,  WA)  was  used  to  calculate  the  angular  velocity  in  real  time  from  the  encoder  readings. 

Simulink  (Mathworks,  Natick  MA)  was  used  to  implement  the  real-time  digital  control  system.  The  torque,  position  and 
velocity  signals  were  sampled  at  4  kHz  (PC6259,  National  instruments,  Austin,  TX),  filtered  using  a  linear  phase  low  pass 
filter  with  the  cutoff  at  the  Nyqulst  frequency  and  converted  into  physical  units.  The  position,  a,  and  the  velocity,  a,  were 
then  multiplied  by  the  user-selected  stiffness  and  damping  coefficients,  k„  and  b^  respectively,  and  then  combined  to  give 
the  target  motor  torque,  t.  This  torque  was  then  applied  to  the  motor  using  the  manufacturer-supplied  torque-voltage 
conversion  factor.  The  servo  motor  controller  (Model:  DPRALTE-020B080,  Advanced  Motion  Controls,  Camarillo,  CA) 
employed  an  Internal  PID  controller  operating  at  20  kHz  ensuring  that  the  torque  applied  to  the  motor  remained  accurate, 
independent  of  load,  speed  and  position. 

Although  arbltraiy  plate  Inertia,  I,  can  be  simulated  using  cyber-physical  systems  (see,  for  example.  Hover  et  al.,  1997; 
Mackowski  and  Williamson,  2011)  the  current  work  Is  restricted  to  cases  in  which  the  effective  inertia  was  the  same  as  the 
physical  Inertia  (/  =  Ip,  where  the  subscript  “p"  here,  and  elsewhere  denotes  a  physical  quantity).  As  described  above,  the 
physical  system  does  not  possess  a  physical  torsional  spring  (kp=0)  and  the  effective  stiffness  is  solely  due  to  the  virtual 
spring  provided  by  the  cyber-physical  system  (k  =  ky).  Lastly,  the  effective  structural  damping  is  provided  by  a  combination 
of  inherent  physical  damping  due  to  friction  in  the  bearings,  and  the  virtual  damping  added  by  the  cyber-physical  model 
(b  =  bp  +  bv).  Ideally  the  physical  damping  should  be  small  compared  with  the  virtual  damping,  and  this  will  be  confirmed  In 
the  following  section.  Lastly,  we  note  that  the  energy  “dissipated”  by  the  virtual  damping  system  represents  energy  that, 
with  appropriate  electronic  circuitry,  could  be  hai'vested  and  put  to  practical  use. 
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Dynamical 
Simulation  Loop 
4  kHz  update  rate 


Fig.  3.  The  cyber-physical  system  control  scheme  composed  of  two  feedback  loops:  (1)  a  current  following  loop  with  a  20  kHz  update  rate  internal  to  the  servo 
controller  and  (2)  a  dynamical  simulation  loop  operating  at  4  kHz  calculates  a  target  torque  value  based  on  the  user-specified  torsional  stiffness  and  damping  gains. 

2.4.  Testing  procedures 

Prior  to  any  fluid-structure  experiments,  an  extensive  series  of  validation  experiments  was  performed  to  assess  the 
behavior  of  the  wind  tunnel  model  and  specifically  (i)  to  characterize  the  internal  (physical)  damping  of  the  motor-bearing 
assembly,  bp,  and  (ii)  to  confirm  that  the  cyber-physical  system  was  capable  of  accurately  simulating  the  linear  structural 
dynamics  of  a  lightly  damped  spring-mass  system.  Once  the  suitability  of  the  model  was  confirmed,  a  wide  variety  of 
measurements  were  taken  to  assess  the  fluid-structure  interactions,  and  the  energy  haiT/esting  performance  of  the  cyber¬ 
physical  flat  plate. 

In  each  trial,  the  value  of  torsional  spring  stiffness,  k,  was  varied,  while  the  other  parameters  (i.e.,  free-stream  velocity, 
L/oo,  equilibrium  angle  of  attack,  Uo,  damping,  b,  and  plate's  mass  moment  of  inertia,  1)  were  all  held  constant.  Thus,  the 
torsional  spring  stiffness  served  as  the  variable  to  modulate  the  structural  frequency  and  the  stability  of  the  system.  In  most 
experiments  the  dimensionless  stiffness  was  initially  set  at  a  high  value  and  then  reduced  incrementally.  A  range  of 
damping  coefficients  and  wind  speed  were  considered  for  the  two  different  plates.  The  values  of  stiffness,  damping,  and 
wind  speed  were  chosen  such  that  the  dimensionless  stiffness  and  damping  spanned  roughly  the  same  range  across  all 
experiments.  These  experiments  enabled  us  to  systematically  observe  the  onset  of  instabilities,  to  characterize  the  general 
behavior  of  the  system  and  to  validate  the  inertial  scaling  proposed.  Note  that,  by  directly  controlling  the  structural  stiffness 
and  damping,  k  and  b,  we  are  able  to  change  the  fluid-structure  interaction  parameters,  k*  and  b*,  while  keeping  the 
Reynolds  number  (and  thus  the  pure  fluid  effects)  constant. 

In  post-processing,  the  aerodynamic  torque  was  determined  from  the  equation  of  motion  (Eq.  (1)) 

Tf(t)  =  Ia(t)  +  bpa(t)  +  TM,  (5) 

where  Tm(t)  =  bya{t)  +  kya(t).  To  evaluate  the  terms  in  the  equation,  the  position  and  torque  signals  were  filtered  by 
employing  a  zero-phase,  sixth  order  Butterworth  low-pass  filter  with  a  cut-off  frequency  set  to  10  times  the  motion 
frequency.  Second-order-accurate  finite  differences  were  used  to  evaluate  the  angular  velocity  and  acceleration. 


3.  Results  and  discussion 

3.1.  Cyber-physical  system  validation  tests 

Several  tests  were  performed  to  determine  the  stiffness  and  damping  characteristics  of  the  system.  First,  the  linearity  of 
the  cyber-physical  torsional  spring,  was  confirmed  for  spring  stiffnesses  ranging  from  0  to  4  Nm/rad.  Next,  the  physical 
damping  associated  with  the  motor  and  shaft  bearings,  bp,  was  determined  by  measuring  the  torque  associated  with 
rotating  the  motor  and  shaft  over  a  range  of  constant  angular  velocities.  The  measured  data  is  shown  in  Fig.  4a,  from  which 
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Fig.  4.  (a)  Damping  model  identification.  The  linear  damping  model  is  plotted  in  blue  for  (d  >  0)  and  red  for  (d  <  0).  The  reduced  linear  model  is  plotted  in 
black,  (b)  Comparison  of  the  experimentally  observed  free-oscillation  behavior  with  the  simulation:  initial  displacement=40\  k  =  2.49  N  m/rad, 
b  ^  0.002  N  m  s/rad,  and  I  ^  7.82  ■  10“^  kg  m^. 

we  see  that  the  friction  is  composed  of  both  Coulomb  and  linear  damping  components 

Tm  +  b(.a)  =  bi{a)  +  b2  Sgn{a),  (6) 

where  r^+b  denotes  the  total  frictional  torque  due  to  motor  and  bearing.  The  coefficients,  b]  and  b2,  were  found  by  fitting 
the  linear  regression  line  to  the  experimental  data.  The  resulting  regression  equations  are  provided  below 

Tm  +  i)(“)  =  2.05  ■  10^^(q;)  +  8.3  ■  lO^^sgn(d),  {a  >  0) 

Tm  +  ijfd)  =  1.62  ■  lO^^fd)  — 8.6  ■  lO^^sgnld),  (d  <  0).  (7) 

Although  one  would  expect  the  damping  model  to  be  symmetric,  the  difference  between  positive  and  negative  rotations 
are  small,  and  within  the  range  of  experimental  uncertainty.  As  will  be  shown  later,  the  aerodynamic  torques  experienced 
during  the  aeroelastic  motions  are  at  least  an  order  of  magnitude  larger  than  the  inherent  physical  damping,  and  for  this 
reason,  it  proved  useful  to  define  a  simplified  damping  model  by  assuming  zero  Coulomb  damping;  this  is  shown  as  a  black 
line  in  Fig.  4a.  In  practice  the  use  of  the  simplified  model  yielded  identical  results  to  the  those  obtained  with  the  full 
damping  description,  and  so  the  simplified  model  was  adopted  throughout. 

The  last  characterization  experiment  performed  was  a  series  of  free  oscillation  (“ring  down”)  experiments,  which  were 
performed  for  different  values  of  plate  inertia,  torsional  spring  stiffness  and  damping  coefficients.  The  response  of  the 
system  from  each  trial  was  compared  with  that  of  a  numerical  simulation  using  the  same  parameter  values.  These 
experiments  were  performed  both  with  and  without  the  plate  mounted  so  that  different  moments  of  inertia  could  be  tested, 
and  to  confirm  that  the  aerodynamic  damping  for  these  tests  was  negligible.  In  general,  the  results  of  these  tests  revealed 
excellent  agreement  between  experiment  and  simulation  (Fig.  4b),  although  for  most  values  of  tested,  the  experiment  is 
observed  to  damp  out  slightly  faster  than  the  simulation  at  the  lowest  amplitudes,  a  discrepancy  believed  to  be  attributable 
to  the  simplified  treatment  of  the  Coulomb  friction  discussed  above.  Nevertheless,  all  the  dynamical  motions  that  we 
present  associated  with  fluid-structure  interactions  have  much  larger  angular  displacements,  and  so  these  low  amplitude 
discrepancies  are  deemed  insignificant  in  practice. 

The  mass  moment  of  inertia  of  the  experimental  apparatus  (including  the  contributions  from  the  plate,  servomotor,  shaft 
coupling,  reaction  torque  sensor  and  encoder)  were  determined  from  the  damping  ratio,  and  the  undamped  natural 
frequency  of  the  system,  co„.  The  damping  ratio  was  obtained  by  calculating  the  rate  at  which  the  damped  oscillation 
amplitude  decreased,  which  is  given  by  the  natural  logarithm  of  the  ratio  of  any  two  successive  amplitudes.  This  technique 
is  commonly  referred  to  as  the  method  of  logarithmic  decrement  (for  example,  Dukkipati,  2004).  The  range  of  damping  ratio 
considered  in  the  present  study  is  approximately  0.002  <.*<0.15.  The  natural  frequency  was  calculated  from  the  simple 
relationship,  coij  =  (On\/\  -  where  =  2;r/T(j.  The  natural  period  of  damped  oscillation,  T^,  was  determined  by  calculating 
the  average  period  between  two  successive  peaks.  Once  o/n  and  k  are  known,  the  mass  moment  of  inertia  can  be  determined 
directly  from  the  expression,  =  y/I^I.  Using  this  method,  the  values  of  mass  moment  of  inertia  were  determined  to  be 
3.45  ■  lO^"*  +  0.01  ■  10“^  kg  m^  (c  =  0.1  m)  and  7.82  ■  lO^"*  +  0.02  ■  lO^"*  kg  m^  (c  =  0.15  m).  The  mass  moment  of  inertia  of 
the  motor-shaft  assembly  (excluding  the  plate)  was  1.98  ■  10“^  +  0.01  ■  10“^  kg  m^. 
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3.2.  Static  aerodynamic  moment 

The  pitching  moment  coefficient  about  mid-chord  was  measured  experimentally  using  two  different  sized  plates  by 
varying  the  angle  of  attack  between  approximately  +90°  at  Re  =  9.9- 10^',  as  illustrated  in  Fig.  5a.  The  value  of  ;r/2, 
predicted  by  the  thin-airfoil  theoiy  is  plotted  with  a  black  solid  line.  The  present  results  compare  well  with  previously- 
reported  data  of  Page  and  Johansen  (1927),  Pelletier  and  Mueller  (2000),  and  Amandolese  et  al.  (2013).  A  careful  and  high- 
resolution  measurement  of  the  pitching  moment  at  fixed  angles  is  shown  in  Fig.  5b  illustrating  a  deviation  from  the  thin- 
airfoil  results  of  Cm  =  ira/2.  At  angles  very  close  to  zero  we  note  that  the  pitching  moment  quickly  adopts  a  nonzero  value  for 
any  angle  of  attack  different  from  zero.  This  is  presumably  due  to  the  fact  that  as  soon  as  the  sharp-edge  plate  deviates  from 
a  =  0,  a  small  separation  bubble  forms,  generating  a  finite  pitching  moment.  As  the  angle  increases.  Cm  increases 
monotonically,  although  the  curvature  changes  from  positive  at  small  angles  to  negative  at  higher  incidence  angles,  with  the 
inflection  point  located  at  a  s!  0.07.  The  significance  of  this  shape,  as  well  as  the  dotted  lines  and  red  circles  shown  in  Fig.  5b 
will  be  discussed  in  context  of  the  dynamic  motion  of  the  plate,  in  the  following  section.  Finally,  while  not  shown  here,  we 
also  confirmed  that  the  slope  of  the  moment,  dCm/da,  is  independent  of  Reynolds  number. 


Fig.  5.  (a)  The  evolution  of  static  moment  coefficient,  C^ia),  as  a  function  of  the  angle  of  attack,  a,  at  Re  =  9.9  •  10'*:  black  circle  (c=0.1  m)  and  black  square 
(c=0.15  m);  thickness-to-chord  ration  1%.  Page  and  Johansen  (1927),  Res*  10“*,  thickness-to-chord  ratio=3%  (red  circle).  Pelletier  and  Mueller  (2000), 
Resi  8  ■  10^,  thickness-to-chord  ratio=1.93%  (blue  circle).  Amandolese  et  al.  (2013),  Re  «i2.3  •  10'',  thickness-to-chord  ratio=4.3%  (green  circle).  Thin-airfoil 
theory,  xa/2  (solid  black  line).  A  best-fit  line  through  the  linear  region  with  a  slope  of  1.53  (dashed  line).  The  results  of  Page  and  Johansen  (1927)  and 
Pelletier  and  Mueller  (2000)  were  adopted  from  Amandolese  et  al.  (2013).  (b)  C„,(a)  versus  a  within  a  small  range  of  a  is  plotted  with  a  blue  curve  (only  the 
first  quadrant  is  shown).  Various  k*a  lines,  including  the  critical  torsional  restoring  torque,  1 .37,  are  plotted  with  dotted-lines.  The  thin  airfoil  theory 

prediction  is  plotted  with  a  solid  black  line.  Approximate  locations  of  the  unstable  (open  red  circle)  and  stable  (solid  red  circles)  fixed  points  are  given  by 
the  intersecting  points  between  the  C^(a)  curve  and  the  /f*a  lines.  The  origin  is  the  unstable  fixed  point.  (For  interpretation  of  the  references  to  color  in  this 
figure  caption,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


3.3.  Flow-induced  oscillations 


We  now  turn  to  the  unsteady  response  of  the  plate  and  describe  its  behavior  for  different  values  of  the  scaled  stiffness, 
damping  and  Reynolds  number.  For  each  experiment,  the  plate  was  initially  positioned  parallel  to  the  direction  of  free 
stream  (a  =  0)  and  the  structural  stiffness  was  set  to  a  large  value  as  to  keep  the  plate  stationaiy.  The  structural  stiffness,  k*, 
was  then  decreased  incrementally,  while  keeping  the  free  stream  velocity  (and  thus,  Reynolds  number),  damping,  and  mass 
moment  of  inertia  constant. 


3.3.1.  Amplitude  response  and  hysteresis 

An  overall  phase  map  of  the  aeroelastic  system  response  as  a  function  of  the  inverse  dimensionless  stiffness  is  illustrated 
in  Fig.  6,  Flere  the  amplitude  of  any  unsteady  motion,  Aa,  is  plotted  against  the  inverse  stiffness,  1  //<*.  The  dynamical 
response  of  the  system  can  be  categorized  broadly  into  four  distinct  regimes:  (1)  for  small  values  of  1  /k*,  the  plate  remains 
stationary  and  nearly  parallel  to  the  flow.  (II)  As  the  stiffness  decreases  below  a  critical  value,  a  region  of  small 
amplitude  limit-cycle  oscillations  (LCOs)  is  obseiA/ed.  These  oscillations  are  asymmetric  (e.g.  Fig.  7,  Region  II),  meaning  that 
the  midpoint  of  the  LCOs  (i.e.,  (max(a)+min(a))/2)  is  non-zero.  (Ill)  As  the  structural  stiffness  is  reduced  further,  the  system 
experiences  a  second  instability  and  the  pitching  amplitude  grows  markedly,  reaching  a  state  characterized  by  pitching 
amplitudes  of  Aasjf.S,  and  symmetric  about  a  =  0.  (IV)  Finally,  as  the  stiffness  is  reduced  even  further,  the  coupling 
between  the  fluid  and  structure  breaks  down  and  the  plate  stops  oscillating,  adopting  a  stable  deflected  state,  perpendicular 
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Inverse  Stiffness,  1/A: 


Fig.  6.  An  overall  phase  map  of  the  aeroelastic  system  response  as  a  function  of  the  inverse  dimensionless  stiffness  for  Re  —  9.9  ■  10“*,  b*  =  0.13,  and 
I*  =  10.6.  The  local  bifurcation  behavior  of  the  system  is  illustrated  in  blue.  Asymmetric  LCOs  are  created  as  the  system  encounters  a  subcritical  pitchfork 
bifurcation  at  the  static  torsional  divergence  condition,  followed  by  a  saddle-node  bifurcation.  Large  symmetric  LCOs  are  generated  as  a  result  of  a 
subcritical  Hopf  bifurcation,  coupled  to  a  saddle-node  bifurcation.  Black  arrows  indicate  the  orientation  of  the  hysteresis.  Three  distinct  dynamical  regimes 
are  labeled  as  1  (steady),  11  (small  asymmetric  LCOs)  and  ill  (large  symmetric  LCOs),  (For  interpretation  of  the  references  to  color  in  this  figure  caption,  the 
reader  is  referred  to  the  web  version  of  this  paper.) 


Fig.  7.  A  magnified  view  of  the  Regions  1  and  11,  illustrating  a  quasi-steady  deflected  state  followed  by  an  unsteady  state  characterized  by  small-amplitude 
asymmetric  LCOs.  A  small  region  of  hysteresis  (0.67  <  1/k*  <  0.73)  is  portrayed.  The  error  bars  denote  the  standard  deviation  of  the  oscillation  amplitude. 


to  the  flow  (this  regime  is  not  portrayed  in  Fig.  6  or  7).  Similar  dynamical  regimes  have  been  reported  by  Dimitriadis  and  Li 
(2009),  who  studied  the  stall  flutter  behaviors  of  a  symmetric  airfoil  with  pitch-plunge  degrees  of  freedom  and  also  by  Kim 
et  al.  (2013)  who  examined  a  self-excited  flapping  behavior  of  the  inverted  flexible  flag  model  with  a  free  leading  edge  and  a 
clamped  trailing  edge.  We  discuss  each  of  these  regimes  in  more  detail  in  the  next  few  paragraphs. 

In  Region  I,  when  1  /k*  is  small,  the  plate  remains  stationary  (Aa  =  0).  As  1  /k*  increases,  the  plate  deflects  slightly  while 
maintaining  a  stable  state  (Fig.  7,  Region  1).  This  subtle  deflection  with  decreasing  structural  stiffness  can  be  understood  in 
conjunction  with  the  static  pitching  moment  (Fig.  5b).  At  any  finite  value  of  k*,  the  plate  adopts  an  angle  such  that  the  static 
aerodynamic  moment  is  balanced  by  the  restoring  torque  of  the  torsional  spring,  k*a.  Since  Cm(a)  is  concave  up,  as  k* 
decreases,  this  stable  fixed  point  moves  to  higher  values  of  a  (Fig.  5b). 

At  a  critical  value  of  1  s!  0.73,  the  plate  enters  Region  11,  in  which  the  plate  exhibits  low-amplitude  asymmetric  LCOs 
(Fig.  7).  At  this  critical  stiffness  value,  based  on  the  static  pitching  characteristics  (Fig.  5b),  the  location  of  the  fixed  point  is 
expected  to  Jump  from  a  Ri  0.03  to  0.11  and  this  agrees  very  well  with  the  experimental  observation.  The  initial  LCO  has  low 
amplitude,  but  it  increases  smoothly  as  1  /k*  continues  to  increase.  The  critical  stiffness  at  which  we  expect  to  see  this 
transition  to  small  amplitude  LCO  can  be  predicted  using  the  classical  quasi-steady  Theodorsen  analysis  of  static  torsional 
divergence  (see,  e.g.,  Bisplinghoff  et  al.,  1957)  in  which  the  dynamical  equation  for  the  motion  of  the  plate  (Eq.  (1))  reduces 
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Fig.  8.  The  evolution  of  angular  position  as  a  function  of  the  dimensionless  time,  Ugot/c,  during  a  transition  from  small  amplitude  asymmetric  LCDs  to  large 
amplitude  symmetric  LCOs.  In  this  case,  i?e  =  9.9’  10“^,  b*  =  0.03.  /*  =  10.6  and  ao  =  0.  The  change  in  vibrational  mode  is  triggered  by  reducing  the 
dimensionless  stiffness,  k*,  from  0.91  to  0.86.  The  vertical  dotted-line  shows  the  instant  at  which  the  structural  stiffness  was  reduced. 


Fig.  9.  Measured  pitching  frequency  versus  the  structural  frequency  (both  normalized  by  the  chord  and  free  stream  velocity)  drawn  from  a  wide  range  of 
experiments;  Uoo  ^  10 25  m/s,  /*  ^4.8  and  10.6,  k*  ^  0.3 1.1,  and  b*  =  0.01 0.21. 

to  =  dCm/da.  The  original  Theodorsen  analysis  (based  on  thin-airfoil  theoiy)  predicts  =  7r/2  =  1.57,  although  if  we 
use  our  experimentally  measured  value  of  Cm(a)  (Fig.  5b),  we  estimate  sj  1.4,  which  is  in  very  good  agreement  with  the 
measured  instability  threshold  of  (Fig.  6).  The  analysis  also  predicts  that  k*^i  is  independent  of  damping  and 

inertia,  and  this  was  confirmed  in  our  experiments  (as  long  as  the  damping  is  not  too  large).  From  a  dynamical  systems 
perspective,  the  system  experiences  a  subcritical  pitchfork  bifurcation  at  the  static  torsional  divergence.  At  this  instance  the 
stable  fixed  point  transforms  into  a  saddle  point  and  two  additional  unstable  equilibria  (spiral  points)  are  created  at  positive 
and  negative  values  of  a.  Subsequently  the  system  undergoes  a  saddle-node  bifurcation  that  ultimately  gives  rise  to  a  stable 
asymmetric  LCO  (Strogatz,  1994).  We  believe  that  the  interaction  between  the  static  divergence  and  dynamic  stall  is 
responsible  for  the  generation  of  these  small  asymmetric  LCOs;  this  assertion  is  reinforced  by  the  fact  that  the  plate 
oscillates  in  and  out  of  the  static  stall  angle  of  a  Ri  0.12.  Our  experimental  observations  and  theoretical  analysis  are  also  in 
agreement  with  those  of  Dimitriadis  and  Li  (2009)  in  which  the  onset  of  small  asymmetric  LCOs  was  also  seen  just  past  the 
static  divergence  condition.  The  pitching  frequency  of  the  small-amplitude  asymmetric  LCOs  (Fig.  7b,  Region  11)  was  difficult 
to  measure  accurately  due  to  the  fact  that  the  oscillations  are  irregular.  However,  to  the  extent  that  they  could  be 
characterized,  they  were  observed  to  be  approximately  20-40%  lower  than  the  natural  structural  frequency  of  the  plate. 

As  1  /k*  continues  to  rise,  we  observe  an  abrupt  transition  to  large-amplitude  symmetric  LCOs  which  we  term  Region  111 
(Fig.  6).  The  transition  occurs  immediately  upon  crossing  the  critical  value  of  1  /k*,  with  the  oscillation  amplitude  growing 
exponentially  in  time  from  its  small-amplitude  state  in  Region  11  to  the  large-amplitude  state  characteristic  of  Region  111 
(Fig.  6).  A  representative  transient  response  is  depicted  in  Fig.  8,  where  the  exponentially  growing  amplitude  is  limited 
presumably  by  unsteady  separated  flows  at  large  angle  of  attack.  This  self-limiting  amplitude  oscillation  is  a  unique 
characteristic  of  the  dynamic  stall  flutter  phenomenon,  which  is  caused  by  the  periodic  separation  and  reattachment  of  the 
flow  around  the  unsteady  plate  immersed  in  a  uniform  free-stream  (McCroskey,  1982).  As  the  phase  portrait  indicates 
(Fig.  6),  this  transition  bears  all  the  hallmarks  of  a  subcritical  Hopf  bifurcation  (Strogatz,  1994),  illustrating  an  abrupt 
transition  to  the  large  amplitude  state,  and  a  hysteretic  return  to  the  stationary  state  if  one  were  to  increase  the  stiffness 
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(i.e.  decrease  1  /k*)  while  in  large  amplitude  regime.  As  shown  in  Fig.  6,  large  amplitude  LCDs  will  persist  until  the  system 
encounters  a  saddle-node  bifurcation  at  1  /k*  0.3,  where  the  stable  and  unstable  LCOs  coalesce.  A  similar  hysteretic 

behavior  has  been  also  obsei'ved  by  a  number  of  investigators  (e.g.,  Razak  et  al.,  2011;  Amandolese  et  al„  2013;  Kim  et  al., 
2013),  although  the  ranges  of  the  hysteresis  reported  in  these  studies  are  generally  smaller  than  what  we  observed  in  the 
present  study.  In  Region  111,  the  oscillations  are  symmetric  about  a  =  0  and  have  a  nearly  sinusoidal  character  with  a  well- 
defined  frequency.  We  find  that,  over  a  very  wide  range  of  velocities,  stiffness  and  (light)  damping,  the  pitching  frequency  is 
very  close  to,  but  always  slightly  lower  than,  the  natural  structural  frequency  of  the  plate  (Fig.  9).  In  this  regime,  any  large- 
scale  vortex  formation  and  shedding  that  must  accompany  these  huge  amplitude  excursions  is  driven  by  the  structural 
dynamics  but  does  not  influence  the  overall  pitching  frequency.  An  explanation  for  the  observation  that  the  pitching 
frequency  is  consistently  below  the  natural  frequency  might  be  that  the  vortex  that  forms  at  the  leading  edge  of  the  plate  as 
it  pitches  up  applies  an  aerodynamic  torque  to  the  plate  and  serves  to  reduce  the  effective  stiffness  of  the  torsional  spring, 
thus  lowering  the  effective  natural  frequency.  This  hypothesis  will  be  revisited  in  Section  3.3.2. 

3.3.2.  Aerodynamic  torque,  energy  harvesting  and  Reynolds  number  effects 

Fig.  10  shows  a  representative  projection  of  the  phase  plane  trajectory,  plotted  over  approximately  10  oscillation  cycles.  A 
typical  phase  plot  of  the  moment,  rp  versus  the  angle,  a,  exhibits  regions  of  both  positive  and  negative  hysteresis.  A  positive 
hysteresis  implies  net  positive  energy  exchange  between  the  fluid  and  the  plate,  in  which  the  energy  is  extracted  from  the 
fluid  and  dissipated  (or  harvested)  by  the  cyber-physical  system.  In  contrast,  negative  hysteresis  implies  that  the  energy  is 
being  transferred  from  the  system  into  the  wake.  The  measured  static  moment  coefficient  curve  is  drawn  in  blue  for 
comparison. 

Several  features  are  worth  discussing.  Starting  at  the  equilibrium  position  of  (a  =0,Tf  =  0),  we  see  that,  as  the  plate 
pitches  up,  the  torque  rises,  but,  curiously,  at  a  slower  rate  than  that  predicted  by  the  static  Cm(a)  curve.  This  suggests  that  a 
dynamically  generated  flow  structure,  possibly  a  vortex  on  tbe  pressure  side  of  the  plate  that  has  survived  from  the  previous 
cycle,  is  resisting  the  positive  rotation  of  the  plate  with  a  small  clockwise  torque.  This  idea  is  supported  by  the  observation 
(not  shown  here)  that,  for  lower  values  of  k*  where  the  pitching  frequency  is  lower  and  unsteady  effects  presumably  weaker, 
the  dynamic  torque-angle  trajectory  falls  closer  to  the  statically  measured  Cm(a)  curve. 

As  the  plate  continues  to  pitch  up,  the  aerodynamic  torque  continues  to  grow  well  beyond  the  static  stall  angle  finally 
reaching  a  maximum  value  that  far  surpasses  its  static  counterpart.  This  is  a  classic  example  of  “dynamic  stall”,  in  which  the 
force  and  torque  on  the  plate  is  strongly  affected  by  the  formation  and  growth  of  a  strong  leading-edge  vortex  (LEV)  that 
develops  on  the  suction  side  of  the  plate  (e.g.,  McCroskey,  1982;  Ohmi  et  al.,  1991 ;  Baik  et  al„  2012;  Ford  and  Babinsky,  2013; 
Granltind  et  al.,  2013),  The  torque  reaches  a  maximum,  after  which  it  begins  a  precipitous  decline,  starting  roughly  at  a  Ri  0.7, 
suggesting  a  stall  condition  most  likely  caused  by  the  separation  of  the  LEV.  A  precise  definition  of  LEV  “separation”  is 
elusive  in  the  literature,  and  the  question  still  remains  to  be  answered  whether  the  LEV  pinch  off  coincides  with  the 
aerodynamic  force  peak,  as  discussed  by  Granlund  et  al.  (2013).  Despite  this  drop  in  torque,  the  plate  continues  to  rotate 
(due  to  its  inertia),  eventually  coming  to  a  rest  at  a  Ri  1.5.  At  this  position,  we  observe  a  second  increase  in  the  aerodynamic 
torque  which  we  hypothesize  is  due  to  the  development  of  a  second  vortex  forming  as  the  shear  layer  shed  from  the  leading 
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Fig.  10.  Aerodynamic  torque  versus  angular  position  phase  portrait  near  the  peak  power  coefficient.  Re  =  9.9  ■10'',  k"  =  0.78.  b*  =  0.13,  /'*  =  10.6,  and 
ao  =  0.  Curved  arrows  indicate  the  direction  of  the  hysteresis.  The  mean  steady  moment  coefficient  curve  is  also  plotted  for  comparison.  The  maximum 
aerodynamic  torque  induced  by  the  pitching  plate  greatly  exceeds  its  static  counterpart;  the  aerodynamic  torque  continues  to  increase  well  beyond  the 
static  stall  angle  of  af«0.12. 
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edge  rolls  up  in  the  vicinity  of  the  plate's  upper  surface.  At  this  extreme  angle  of  rotation,  the  structure's  torsional  spring 
exerts  a  strong  torque  which  pulls  the  plate  back  towards  its  equilibrium  position.  However,  as  the  plate  starts  its  return,  a 
small  region  of  negative  hysteresis  is  observed  owing  to  the  fact  that  the  plate's  return  motion  is  hampered  by  the 
aerodynamic  counter-torque  induced  by  the  second  vortex.  As  the  angle  of  incidence  becomes  small,  we  anticipate  that  the 
flow  over  the  plate  to  temporarily  reattach  to  the  surface  and  the  aerodynamic  torque  returns  to  zero.  Due  to  the  symmetiy 
of  the  configuration,  the  pattern  repeats  during  the  second  half  of  the  cycle. 

The  net  aerodynamic  work  per  oscillation  cycle  can  be  determined  by  calculating  the  area  enclosed  by  the  zf-a  hysteresis 
loop 

I  l‘a{to-\-NJ)  I  1‘to+NT 

'"  =  nL  = 

where  T  is  the  cycle  period,  and  N  is  an  integral  number  of  cycles.  The  integration  was  performed  by  employing  a  trapezoidal 
integration  scheme  with  a  time  step  of  0.25  ms  using  30  s  of  recorded  data  (N  r;  0(10^)).  The  power  coefficient,  Cp,  is  simply 
the  energy  per  cycle,  W,  multiplied  by  the  pitching  frequency  and  appropriately  normalized 

Cp  =  WT/^pUich.  (9) 

Fig.  11  shows  contour  plots  of  the  power  coefficient  at  three  different  Reynolds  numbers.  The  power  coefficient  is 
characterized  over  a  range  of  dimensionless  stiffness  and  damping.  First,  we  note  that  the  magnitude  of  the  power 
coefficient  remains  approximately  independent  of  Reynolds  number,  reaching  a  maximum  value  of  approximately 
1.2%.  Secondly  the  range  of  k*  over  which  power  generation  is  achieved  is  also  approximately  constant  with  Reynolds 
number  with  the  optimal  stiffness  for  power  generation  located  slightly  lower  than  unity.  Both  of  these  results  confirm 
that  the  inertial  scaling  chosen  to  describe  the  physics  of  these  large  scale  oscillations  is  appropriate  for  this  problem. 
Nevertheless,  there  is  an  observable  effect  of  Reynolds  number.  As  Re  increases,  larger  values  of  structural  damping,  b*, 
can  be  employed  while  still  maintaining  large  scale  motion  to  extract  energy  from  the  flow.  It  is,  nevertheless 
interesting  that  despite  using  larger  values  of  b*,  the  power  coefficient  does  not  increase  appreciably.  This  is  likely  due 
to  the  fact  that  although  more  energy  per  cycle  is  being  harvested  (recall  that  the  structural  damping  is  the  mechanism 
by  which  energy  is  harvested  from  the  fluid),  the  pitching  frequency  and  the  amplitude  are  reduced.  The  phase 
portraits  of  torque,  zf,  versus  position,  a,  associated  with  two  different  Reynolds  numbers  and  two  different  damping 
coefficients  are  shown  in  Fig.  12.  In  each  case,  k*  and  I*  were  held  constant.  The  corresponding  cases  are  shown  using  a 
green  triangle  and  circle  in  Fig.  11.  For  a  moderately  low  damping  case  (Fig.  12a,  and  the  circles  in  Fig.  11),  one  can  see 
that  the  shape  and  the  size  of  the  hysteresis  trajectories  are  very  similar;  consequently  the  dimensionless  pitching 
frequency,  amplitude,  and  power  coefficient  scale  very  well  at  two  different  Re  values.  On  the  contrary,  the  Reynolds 
number  effect  appears  to  become  more  prominent  at  a  higher  dimensionless  damping  value,  as  depicted  in  Fig.  12b 
(the  triangles  in  Fig.  11).  We  can  clearly  see  that  not  only  the  maximum  pitching  amplitude  is  increased,  but  also  the 
plate  experiences  slightly  larger  aerodynamic  torque.  The  high  Reynolds  number  case  generates  a  much  higher  power 
coefficient  (0.0092  vs.  0.0054)  and  has  a  much  higher  pitching  frequency  (7.1  Hz  vs.  4.5  Hz).  Since  the  structural 
characteristics  of  the  system  are  identical,  the  difference  between  the  two  system's  behavior  must  lie  in  the  fluid 
behavior  and  faster  generation  of  a  stronger  and  more  stable  LEV  is  likely  responsible  for  the  augmentation  in  the 
aerodynamic  torque  and  pitching  amplitude.  Velocity  measurements  are  needed  to  properly  identify  these  effects, 
something  beyond  the  scope  of  this  paper. 


Damping,  b’  Damping,  b~  Damping,  b’ 

Fig.  11.  The  power  coefficient  contour  plots  for  various  Reynolds  numbers:  (left)  ^6  =  9.910“^,  (center)  Re  =  1.110^,  and  (right)  Re=1.310^.  The 
normalized  mass  moment  of  inertia,  I*,  is  10.6.  The  o  and  represent  low  and  high  damping  cases,  respectively,  whose  corresponding  aerodynamic  torque 
versus  angular  position  phase  plane  trajectories  are  illustrated  in  Fig.  12.  The  contours  are  constructed  using  linear  interpolation. 
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Fig.  12.  Changes  in  the  normalized  aerodynamic  torque  versus  angular  position  hysteresis  curves  at  two  distinct  Reynolds  numbers  (Re  =  9.9  ■  1  O'*  and 
Re  =  1.3  ■  10^)  for  two  different  values  of  the  dimensionless  damping:  (a)  b*  =  0.11  and  (b)  b*  =  0.23.  In  both  cases,  the  dimensionless  mass  moment  of 
inertia  and  torsional  stiffness  were,  respectively,  1*  =  10.6  and  k*  =  0.68.  The  data  are  plotted  over  several  cycles.  Cases  (a)  and  (b)  are  designated  with  o 
and  respectively,  in  Fig.  11. 


3.3.3.  Energy  harvesting  enhancement  using  a  nonlinear  (Cubic)  spring 

Perhaps  the  most  appealing  aspect  of  the  cyber-physical  system  is  that  it  enables  the  user  to  easily  explore  the  impact  of 
structural  nonlinearities  (i.e,  nonlinear  elastic  and  damping  forces)  on  the  fluid-structure  interactions.  In  this  section,  we 
demonstrate  the  implementation  of  a  simple  nonlinear  cubic  spring,  and  show  how  it  can  be  utilized  to  improve  the  energy 
extraction  capacity  of  the  system.  The  nonlinear  elastic  torque,  is  modeled  as 

Tv  =  (fO) 

where  )ci  and  K3  are  the  linear  and  cubic  spring  constants,  respectively.  In  this  example,  we  wish  to  model  a  hardening  cubic 
spring  with  values  of  kx  =  0.58  and  K3  =  0.42.  In  this  way,  the  nonlinear  spring  has  similar  characteristics  to  the  original 
linear  spring  for  small  angular  displacement  (Fig.  13a).  With  the  installment  of  the  cubic  spring,  the  dimensionless  pitching 
frequency  and  the  power  coefficient  increased  by  37.5%  and  57.0%,  respectively.  The  phase  portraits  of  t/ versus  a  associated 
with  the  two  distinct  springs  are  presented  in  Fig.  13b  and  can  be  used  to  explain  these  changes.  One  can  clearly  observe 
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Fig.  13.  (a)  Dimensionless  elastic  torque  as  a  function  of  angular  displacement  for  a  linear  and  nonlinear  springs.  Elastic  torque  was  measured 
experimentally,  (b)  Aerodynamic  torque  versus  angular  position  phase  portraits  for  linear  and  cubic  springs.  Re  =  9.9  •  10"^,  b*  =  0.13,  I*  =  10.6,  and  ao  =  0. 
The  values  of  spring  constants  associated  with  a  cubic  spring  (as  defined  by  Eq.  (10))  are  xrs  =  0.42  and  xri  =0.58.  The  linear  spring  constant  equals  kj. 
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that,  by  implementing  a  simple  cubic  spring,  the  regions  of  negative  hysteresis  typically  observed  during  the  pitch  reversal 
have  been  significantly  reduced,  and  the  maximum  pitching  amplitude  has  been  decreased  by  roughly  12%.  More 
importantly,  by  increasing  the  elastic  restoring  force  at  maximum  rotation,  the  plate  is  forced  to  spring  back  faster  (hence 
the  higher  pitching  frequency)  and  we  hypothesize  that  this  reduces  the  time  available  for  the  second  (and  detrimental) 
vortex  to  form  at  the  leading-edge,  thus  eliminating  the  negative  hysteresis.  This  hypothesis  is  yet  to  be  confirmed  (using 
velocity  measurements). 


4.  Concluding  remarks 

The  work  presented  here  demonstrates  a  number  of  advances.  We  have  successfully  demonstrated  the  use  of  a  cyber¬ 
physical  system  to  accurately  and  easily  simulate  a  wide  range  of  structural  dynamics.  In  the  current  work  we  have  not 
implemented  any  virtual  modification  of  the  plate's  inertia,  and  have  restricted  the  study  to  a  single  degree  of  freedom 
motion.  However,  both  these  are  easy  to  implement,  and  despite  the  initial  complexity  of  the  cyber-physical  system  setup, 
the  key  benefit— that  of  rapid  and  extensive  access  to  the  entire  fluid-structural  parameter  space— is  the  abiiity  to  perform  a 
very  deep  study  of  the  physics  of  fluid-structure  interactions.  Unlike  previous  cyber-physical  systems  (Hover  et  al.,  1997;  Lee 
et  al.,  2011;  Mackowski  and  Williamson,  2013)  that  operate  in  water  at  low  frequencies,  our  system  operates  in  air  at  high 
operational  frequencies  and  with  much  larger  and  more  violent  structural  motion. 

Nevertheless,  the  purpose  of  implementing  the  cyber-physical  system  is  not  only  as  a  technical  demonstration,  but  rather 
to  enable  the  exploration  and  detailed  characterization  of  fluid-structure  interactions  and  here  too  we  have  presented 
several  novel  insights.  Although  the  onset  of  the  small-amplitude  nonlinearly  stabilized  oscillations  (Region  II)  has  been 
widely  explored,  dating  back  to  Theodorsen,  the  large-amplitude  motion  (Region  III)  has  not  been  well  characterized,  and 
the  current  data  provides  valuable  insight  into  this  regime.  The  pitching  moment  associated  with  the  large  scale  unsteady 
motion  is  greatly  enhanced,  with  torques  nearly  three  times  as  large  as  those  associated  with  quasi-steady  motion.  These 
forces  are  achieved,  presumably,  by  the  dynamics  of  the  leading-edge  vortex  that  forms  at  the  sharp  leading  edge  of  the 
plate,  and  the  cyber-physical  system  can  be  tailored  using  the  characteristics  of  the  torsional  spring,  to  modify  the 
interaction  between  the  vortex  and  the  structure.  Furthermore,  by  monitoring  the  phase  plane  representation  of  the 
structural  motion,  we  can  follow  the  formation  time  and  separation  dynamics  of  these  leading-edge  structures  and  identify 
timescales  as  a  function  of  Reynolds  number.  Finally,  from  the  perspective  of  energy  harvesting,  with  efficiencies  of  only  one 
or  two  percent,  the  pitching  plate  is  a  poor  choice  if  one  seeks  efficient  extraction  of  energy  from  the  free  stream.  The 
pitching-plunging  plate  studied  by  several  other  researchers  (see,  e.g.,  Dumas  and  Kinsey,  2006;  Simpson  et  al.,  2008;  Peng 
and  Zhu,  2009)  is  much  more  effective  since  it  can  take  advantage  of  the  aerodynamic  lift  during  the  heave  motion  to  extract 
significant  energy  from  the  flow.  It  is  not  our  intent  to  propose  an  efficient  energy  harvesting  system,  but  rather  to  gain 
perspective  on  the  dynamics  of  energy  flow  in  a  simplified  physical  system,  and  here  a  key  insight  gained  from  the  present 
study  is  the  role  that  the  torsional  spring  stiffness  plays  in  establishing  feedback  between  the  fluid  and  the  structure,  the 
role  of  damping  in  moderating  that  feedback,  and  the  appropriate  scaling  of  the  stiffness  and  damping  necessary  for 
maximum  energy  extraction  performance  for  a  given  flow  condition.  A  hysteretic  behavior  of  an  aeroelastic  energy 
harvesting  system  is  beneficial  in  real  applications  as  it  serves  to  expand  the  operational  regime  for  power  generation.  The 
implementation  of  such  devices  in  real  applications  poses  enormous  challenges  owing  to  difficulty  in  tuning  the  system’s 
structural  characteristics  to  realize  stable  and  efficient  power  generation.  One  possible  solution  is  to  implement  a  cyber¬ 
physical  system,  coupled  with  an  optimization  scheme,  as  a  means  to  automate  the  modulation  of  the  structural  parameters. 

Moving  forward,  the  insight  gained  from  the  phase  portraits,  in  particular  with  regard  to  the  behavior  of  the  separated 
flow  structures,  must  be  confirmed  and  quantified  directly,  and  the  next  phase  of  the  research  will  include  velocimetiy 
(using  Particle  Image  Velocimetry),  coupled  to  measurements  of  the  plate  motion.  We  believe  that  the  results  from  these 
experiments  will  shed  more  insight  on  the  physics  of  the  fluid-plate  interactions  and  the  process  of  energy  transfer  from  the 
fluid  to  the  structure. 
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We  report  on  the  dynamics  of  the  formation  and  growth  of  the  leading-edge  vortex  and 
the  corresponding  unsteady  aerodynamic  torque  induced  by  large-scale  flow-induced  os¬ 
cillations  of  an  elastically  mounted  flat  plate.  All  experiments  are  performed  using  a 
high-bandwidth  cyber-physical  system,  which  enables  the  user  to  access  a  wide  range 
of  structural  dynamics  using  a  feedback  control  system.  A  series  of  2D  particle  image 
velocimetry  is  carried  out  to  characterize  the  behaviour  of  the  separated  flow  structures 
and  its  relation  to  the  plate  kinematics  and  unsteady  aerodynamic  torque  generation. 
By  modulating  the  structural  properties  of  the  cyber-physical  system  we  systematically 
analyse  the  formation,  strength,  and  separation  of  the  leading-edge  vortex,  and  the  de¬ 
pendence  on  kinematic  parameters.  We  demonstrate  that  the  leading-edge  vortex  growth 
and  strength  scale  with  the  characteristic  shear-layer  feeding  velocity  and  that  a  poten¬ 
tial  flow  model  using  the  measured  vortex  circulation  and  position  can,  when  coupled 
with  the  steady  moment  of  the  flat  plate,  accurately  predict  the  net  aerodynamic  torque 
on  the  plate.  Connections  to  previous  results  on  optimal  vortex  formation  time  are  also 
discussed. 

Key  words:  Leading-edge  vortex,  flow-induced  oscillations,  cyber-physical  system 


1.  Introduction 

Characterization  of  unsteady  leading-edge  vortex  (LEV)  dynamics  associated  with 
separation  from  bodies  at  high  angles  of  attack  has  been  a  problem  of  scientific  interest 
for  decades.  The  importance  of  LEV  formation  (the  hallmark  of  dynamic  stall)  in  relation 
to  force  augmentation  is  well  documented  in  the  archival  literature,  often  in  the  context 
of  biological  flight  (Ellington  et  al.  1996),  the  dynamics  of  highly  maneuverable  aircraft 
(01  et  al.  2010)  and,  more  recently,  energy  conversion  applications  (Peng  &  Zhu  2009; 
Ramesh  et  al.  2015;  Onoue  et  al.  2015). 

A  number  of  studies  that  imposed  predefined  wing  kinematics  to  characterize  the  LEV 
formation  have  demonstrated  that  the  behaviour  of  the  separated  flow  structures  is  quite 
sensitive  to  variations  in  wing  kinematics  (Kim  &  Gharib  2010;  Baik  et  al.  2012;  Ozen 
&  Rockwell  2012),  pivot  axis  location  (Yu  &  Bernal  2013;  Granlund  et  al.  2013),  wing 
geometry  (Yilmaz  &  Rockwell  2012;  Hartloper  &  Rival  2013),  and  Reynolds  numbers 
(Jantzen  et  al.  2014).  Eor  this  reason,  it  poses  an  enormous  challenge  to  find  unifying 
scaling  principles  for  the  LEV  formation  time,  strength  and  stability  that  can  be  gen¬ 
eralized  to  arbitrary  problems.  One  promising  approach  uses  the  so-called  “universal 
vortex  formation  number” (originally  introduced  by  Gharib  et  al.  (1998)),  which  can  be 
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extended  to  a  wide  variety  of  vortical  flows;  Dabiri  (2009)  defined  the  formation  number 
based  on  the  feeding  shear-layer  velocity  and  the  characteristic  length  of  the  vortex  gen¬ 
erator.  The  premise  of  this  universal  vortex  formation  number  is  that  the  dimensionless 
dynamic  formation  time  corresponding  to  the  generation  of  a  maximum-strength  vor¬ 
tex  remains  constant  at  a  value  of  about  4,  at  which  time  the  vortex  rejects  additional 
entrainment  of  vorticity  from  the  shear  layer  and  pinches  off.  Measurements  by  Rival 
et  al.  (2009)  have  demonstrated  that  the  formulation  of  Dabiri  (2009)  provides  a  good 
scaling  for  the  maximum  LEV  circulation  generated  by  an  aerofoil  executing  sinusoidal 
plunging  motions.  The  use  of  feeding-shear  layer  velocity  to  scale  LEV  circulation  seems 
to  be  a  logical  choice  given  that  the  vortex  growth  is  dictated  by  the  transport  of  circu¬ 
lation  flux  into  the  vortex  core  through  the  feeding  shear-layer,  as  suggested  by  Sattari 
et  al.  (2012).  More  recently,  Kriegseis  et  al.  (2013)  measured  the  LEV  circulation  and 
unsteady  forces  acting  on  a  flat  plate  undergoing  accelerating  plunging  and  towing  mo¬ 
tions  at  a  fixed  angle  of  attack.  The  authors  have  demonstrated  that  the  LEV  circulation 
and  hydrodynamic  forces  can  be  effectively  normalized  by  the  characteristic  shear-layer 
velocity,  Ueff  =  Uoo  +  h,  where  h  designates  the  towing  or  plunging  velocity.  Buchholz 
et  al.  (2011)  has  proposed  an  alternative  scaling  for  the  total  spanwise  circulation  shed 
by  a  periodically  pitching  rectangular  plate  about  its  leading  edge.  This  scaling  parame¬ 
ter  is  formulated  based  on  a  dynamic  pressure  associated  with  the  maximum  transverse 
velocity  of  the  trailing  edge  and  the  plate  geometry.  The  authors  have  demonstrated  that 
their  scaling  parameter  can  be  used  to  satisfactorily  collapse  the  vortex  circulation  data 
for  different  aspect  ratios  and  pitching  amplitudes.  However,  it  still  remains  to  be  veri¬ 
fied  whether  the  aforementioned  scaling  relations  are  robust  enough  to  collapse  the  LEV 
circulation  histories,  both  in  magnitude  and  time,  throughout  the  entire  vortex  growth 
phase. 

While  these  studies  have  used  proscribed  kinematics  to  generate  LEVs,  in  real  sys¬ 
tems,  the  vortex  structures  are  often  a  result  of  fluid  structure  interactions,  in  which  the 
structural  dynamics  of  an  elastically-supported  aerodynamic  surface  (for  example  a  flat 
plate)  are  driven  by  the  growth,  evolution  and  shedding  of  strong  LEVs.  As  discussed 
in  Onoue  et  al.  (2015),  such  aeroelastic  interactions  often  result  in  non-destructive  limit- 
cycle  oscillations  (LCOs)  with  a  nearly  constant  amplitude  and  frequency,  which  can 
potentially  serve  as  the  basis  for  the  passive  power  generation.  While  previous  studies  on 
dynamic-stall  flutter  have  provided  a  wealth  of  information  with  regard  to  the  instabil¬ 
ity  boundaries,  local  bifurcation  behaviours,  and  amplitude  of  LCOs  with  variations  in 
flow  speed  and  structural  characteristics  (Amandolese  et  al.  2013;  Poirel  &  Mendes  2014; 
Onoue  et  al.  2015),  a  direct  connection  between  the  large  amplitude  oscillations  and  the 
vortical  structures  generated  by  the  violent  motion  of  the  structure  has  not  yet  been 
made.  Onoue  et  al.  (2015)  measured  the  aerodynamic  torque  as  a  function  of  the  plate 
motion,  from  which  a  plausible  description  of  vortex  formation,  growth  and  separation 
can  be  inferred,  but  without  measurements  of  the  velocity  field,  it  is  impossible  to  develop 
robust  aerodynamic  models  that  allow  us  to  characterize,  and  thus  predict,  the  detailed 
behaviour  of  large-scale  fluid-structure  interactions  (and  positive  energy  harvesting). 

The  measurements  and  analyses  presented  in  the  current  manuscript  provide  this  valu¬ 
able  connection  between  the  vortex  dynamics  and  aerodynamic  loads,  and  in  this  study, 
we  take  advantage  of  the  cyber-physical  experimental  setup  in  conjunction  with  syn¬ 
chronized  PIV  to  carry  out  a  detailed  study  of  the  flow  physics  associated  with  a  rapidly 
pitching  plate  in  a  uniform  airflow.  As  described  in  Onoue  et  al.  (2015),  the  cyber-physical 
system  is  one  in  which  the  structural  dynamics  (torsional  spring  stiffness,  damping  and 
mass  moment  of  inertia)  are  emulated  using  a  digital  controller  rather  than  by  physi¬ 
cal  hardware.  By  being  able  to  directly  control  the  structural  properties  of  the  system. 
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Figure  1.  (a)  schematic  of  an  oscillating  flat  plate  with  a  pitch  axis  located  at  the  mid  chord. 
The  dashed  horizontal  line  designates  the  equilibrium  angle  of  attack,  ao=0.  (b)  A  computer- 
controlled  servo  motor  rotates  the  shaft-plate  assembly  by  applying  a  specific  torque  calculated 
via  numerical  simulation  to  achieve  the  desired  structural  dynamics.  See  Onoue  et  al.  (2015)  for 
detailed  descriptions  of  the  cyber-physical  system  control  algorithm. 

the  user  can  systematically  tune  the  fluid-structure  interaction  parameters  while  com¬ 
pletely  isolating  the  effects  of  Reynolds  number,  and  hence  the  pure  fluid  effects.  The 
use  of  cyber-physical  system  was  first  demonstrated  by  Hover  et  al.  (1997),  who  success¬ 
fully  implemented  a  force-feedback  algorithm  to  impose  arbitrary  structural  forces  in  the 
study  of  vortex-induced  vibration  (VIV)  of  marine  cables.  In  recent  years,  several  cyber¬ 
physical  systems  have  been  successfully  employed  in  VIV  energy  harvesting  investigations 
(Lee  &  Bernitsas  2011;  Mackowski  &  Williamson  2013).  While  all  of  these  cyber-physical 
systems  have  been  designed  to  be  operated  in  water  at  relatively  low  frequencies,  our 
cyber-physical  system  operates  in  air  (where  the  plate  inertia  is  more  important),  at 
higher  frequencies  and  with  much  more  violent  structural  motions  (Onoue  et  al.  2015). 

The  primary  objectives  of  the  current  study  are  two-fold:  (1)  to  determine  the  scaling 
relation  for  the  LEV  growth  rate  and  circulation  on  a  rapidly  pitching  flat  plate,  and 
(2)  to  elucidate  the  linkage  between  the  LEV  dynamics  and  the  aerodynamic  forces  and 
moment  that  the  flow  imparts  on  the  elastic  structure. 


2.  Experimental  techniques 

Experiments  were  conducted  in  a  closed-return  wind  tunnel  with  a  test  section  of  0.6 
m  in  width,  0.6  m  in  depth,  and  1.2  m  in  length.  The  turbulence  level  in  the  wind  tunnel 
was  previously  measured  to  be  less  than  0.1%  of  the  free-stream  velocity  (Song  et  al. 
2008).  Figure  1(a)  illustrates  a  flat  plate,  with  chord,  c,  span,  h  and  thickness,  6,  free  to 
rotate  about  its  mid-chord  in  a  uniform  airflow  of  magnitude,  Uoo,  and  density,  p/.  In 
the  present  study,  the  Reynolds  numbers  based  on  the  chord  length.  Re  =  Uooc/v,  were 
in  the  order  of  O{10'^)  to  O{10^),  where  u  is  the  kinematic  viscosity  of  the  fluid.  Since 
the  chord  represented  no  more  than  17%  of  the  test  section  width,  blockage  and  side- wall 
boundary  layer  effects  were  deemed  insignificant. 

Figure  1(b)  shows  a  schematic  of  the  experimental  configuration  used  in  the  current 
study,  illustrating  the  flat  plate  mounted  to  a  virtual  mass-spring-damper  system  in  the 
wind  tunnel  test  section.  Instead  of  using  a  physical  torsional  spring  and  damper,  the 
structural  stiffness  and  damping  characteristics  of  the  pitching  plate  were  controlled  using 
a  cyber-physical  system  that  consists  of  a  computer-controlled  servomotor  (SM233AE, 
Parker  Compumotor,  Rohnert  Park,  CA),  torque  sensor  (TFF400,  Futek,  Irvine,  CA), 
and  rotary  encoder  (HB6M  Optical  Encoder,  US  Digital,  Vancouver,  WA)  coupled  to  a 
real-time  digital  controller  operating  in  a  feedback  loop.  The  cyber-physical  system  also 
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sCMOS  camera  system 


Figure  2.  Experimental  setup  for  the  PIV  measurements  in  the  streamwise  plane  at  the 

midspan  location. 


allows  for  the  virtual  modification  of  the  structural  mass  moment  of  inertia.  The  digital 
control  system  (Simulink,  Mathworks,  Natick  MA)  samples  the  angular  position  signal, 
a(f),  from  the  optical  encoder  at  4  kHz,  and  numerically  calculates  the  velocity,  a{t)  and 
acceleration,  a{t),  in  real-time  via  custom-designed  FIR  differentiators  with  a  linear- 
phase  response  in  the  pass-band.  These  signals  are  then  multiplied  by  the  user-specified 
stiffness,  fc„,  damping,  by,  and  mass  moment  of  inertia,  ly,  coefficients,  respectively,  where 
the  subscript  v  denotes  a  virtual  quantity.  These  signals  are  combined  to  determine  the 
target  motor  torque,  T^it)  =  Iya{t)  +  bya{t)  +  kya{t),  which  is  applied  to  the  servomotor. 
This  torque  is  measured  using  the  reaction  torque  sensor  with  a  maximum  measurement 
capacity  of  7  N-m.  In  post-processing,  the  aerodynamic  torque  was  determined  from  the 
equation  of  motion:  Tf(t)  =  Ipuit)  -\-  bpa(t)  -|-  where  the  subscript  p  designates  a 

physical  quantity.  A  series  of  free-oscillation  tests  were  conducted  in  air  to  deduce  the 
value  of  physical  mass  moment  of  inertia.  Ip,  using  the  method  of  logarithmic  decrement. 
These  free-oscillation  tests  exhibit  exponentially-decaying  amplitude  such  that  the  effect 
of  Coulomb  damping  (or  dry  friction)  can  be  neglected.  Therefore,  the  physical  damping 
can  be  simply  modelled  as  a  linear  function  of  velocity,  bpa(t).  The  physical  damping 
associated  with  the  motor  and  bearing  was  determined  by  measuring  the  torque  resulting 
from  rotating  the  motor  and  shaft  over  a  range  of  constant  angular  velocities,  and  fitting 
the  linear  regression  line  to  the  data.  Finally,  to  evaluate  the  above  equation  for  Tf{t), 
the  position  and  torque  signals  were  filtered  by  employing  a  zero-phase,  sixth  order  But- 
terworth  low-pass  filter  with  a  cut-off  frequency  set  to  ten  times  the  pitching  frequency. 
Second-order  accurate  finite  differences  were  used  to  evaluate  the  angular  velocity  and 
acceleration  during  post-processing.  Hereinafter  for  simplicity,  the  effective  torsional  stiff¬ 
ness,  damping  and  mass  moment  of  inertia  are  denoted  without  any  subscript,  as  k,  b, 
and  I,  respectively. 

The  setup  for  the  streamwise  2D  PIV  experiment  is  illustrated  in  figure  2.  All  PIV 
measurements  were  acquired  in  a  streamwise  plane  located  at  the  mid-span  of  the  plate. 
The  laser  sheet  was  generated  using  a  double-pulsed  200  mJ  Nd:YAG  laser  (A=532 
nm)  (EverGreen,  Quantel  USA,  MT)  with  LaVision  sheet  optics.  The  thickness  of  the 
laser  sheet  was  approximately  3  mm.  The  flow  was  seeded  with  oil  droplets  nominally 
1  ^m  in  diameter.  The  mid-span  section  of  the  plate  is  made  out  of  clear  anti-reflective 
acrylic  which  enabled  PIV  measurements  from  both  sides  of  the  structure.  The  particle 
images  were  acquired  with  four  synchronized  sGMOS  cameras  (LaVision)  equipped  with 
50  mm  lenses  (Nikon).  Each  camera  has  a  resolution  of  2560  pixels  x  2160  pixels  and 
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the  frame  rate  for  image  acquisition  was  15  Hz.  Velocity  fields  were  calculated  using 
DaVis  8.6  software  (LaVision)  by  applying  a  sequential  spatial  cross-correlation  with 
multi-pass  iterations  with  decreasing  window  size  (128x128,  2  iterations,  50%  overlap 
to  32x32,  2  iterations,  50%  overlap).  Velocity  helds  acquired  by  four  different  cameras 
were  stitched  together  by  averaging  the  overlapped  regions,  yielding  a  composite  field- 
of-view  of  approximately  0.25  m  by  0.3  m  containing  over  63,000  vectors.  For  each  PIV 
experiment,  approximately  2500  composite  image  pairs  were  acquired,  and  subsequently 
allocated  into  72  equally-spaced  phase  bins  equally  spaced  over  the  oscillation  cycle 
{A(j)=5°).  The  instantaneous  phase  was  calculated  by  taking  the  Hilbert  transform  of 
the  angular  position  time  history  (Khalak  &  Williamson  1999). 

In  order  to  identify  the  vortical  structures  in  the  flow  held,  the  vortex  identihcation 
technique  based  on  the  vortex  swirling  strength  was  employed  (Adrian  et  al.  (2000)).  The 
threshold  was  set  to  10%  of  the  maximum  swirling  vortex  strength  to  isolate  the  vortices 
of  interest  from  the  low-level  background  noise.  This  threshold  was  appropriate  to  exclude 
the  vorticity  conhned  in  the  feeding  shear-layer  in  the  calculation  of  the  primary  vortex 
strength.  The  spanwise  circulation  of  the  primary  vortex  was  calculated  by  integrating 
the  simply  connected  vorticity  patch  about  the  vortex  centroid  using  the  Stokes’  theorem. 
The  results  were  insensitive  to  the  threshold  values  in  the  range  of  0.1±0.05.  Another 
local  vortex  identification  method  (Q-criterion)  was  also  tested,  and  was  found  to  yield 
similar  results.  In  addition,  the  Ti  criterion  (Graftieaux  et  al.  2001)  was  employed  to 
identify  the  location  of  the  LEV  centroid.  In  the  present  study,  a  threshold  of  Ti^O.O 
was  selected,  and  the  vortex  centre  is  found  where  the  maximum  spatial  value  of  Ti 
is  achieved.  This  vortex  identihcation  method  has  been  successfully  implemented  by 
numerous  investigators  (e.g.,  Baik  et  al.  (2012);  Jantzen  et  al.  (2014)). 


3.  Results  and  discussion 

3.1.  Phase-averaged  swirl  fields  and  Tf  vs.  a  phase  portrait 

As  mentioned  earlier,  the  plate  exhibits  very  large  amplitude  limit  cycle  oscillations 
(LCO)  over  a  wide  range  of  structural  parameters,  presumably,  driven  by  the  devel¬ 
opment,  growth  and  separation  of  leading  edge  vortices.  A  representative  sequence  of 
phase-averaged  swirl  helds  at  several  stages  during  the  half  cycle  are  presented  in  hg- 
ure  3,  along  with  the  associated  torque  versus  angular  position  phase  portrait.  As  the 
plate  begins  to  rise  from  its  equilibrium  position,  oq  =  0  (hgure  3:  (1)),  the  aerodynamic 
torque  increases,  but  interestingly  the  dynamic  torque  is  lower  than  that  predicted  by 
the  steady  moment  coefficient  curve  (also  shown  in  hgure  3),  as  discussed  by  Onoue 
et  al.  (2015).  This  is  presumably  due  to  the  apparent  generation  of  the  vortices  at  the 
half-chord  location  as  the  boundary  layer  on  the  suction  side  of  the  plate  separates  from 
the  rotating  shaft  (hgure  3:  (1)).  This  in  turn  generates  a  low  pressure  region  over  the 
rear-half  of  the  plate,  causing  the  centre  of  pressure  to  move  further  aft  of  the  leading 
edge.  Measurements  over  a  wide  range  of  conditions  indicate  that  this  effect  is  weaker 
at  lower  pitching  frequencies,  for  which  cases  the  aerodynamic  torque-angle  trajectory 
falls  closer  to  the  statically-measured  Cm  (a)  curve.  The  torque  offset  at  a=0  can  be 
attributed  to  the  separated  how  on  the  lower  side  of  the  plate  at  the  end  of  the  previous 
cycle,  which  creates  a  pressure  gradient  in  favour  of  generating  an  aerodynamic  torque 
in  the  direction  of  the  wing  rotation. 

As  the  angle  of  attack  increases,  the  aerodynamic  torque  continues  to  grow  well  beyond 
the  static  stall  angle  hnally  reaching  a  maximum  value  nearly  2.5  times  greater  than 
its  static  counterpart  (hgure  3:  (l)-(3)).  This  is  a  classic  example  of  “dynamic  stall” 
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Figure  3.  Snapshots  of  the  phase-averaged  swirl  fields  over  a  half  cycle  and  the  corresponding 
aerodynamic  torque  versus  angular  position  phase  portrait.  Curved  arrows  indicate  the  orienta¬ 
tion  of  the  hysteresis  trajectory.  Regions  of  positive  and  negative  hysteresis  are  shown;  a  positive 
hysteresis  implies  net  positive  energy  exchange  between  the  fluid  and  the  plate  and  a  negative 
hysteresis  describes  the  amount  of  energy  transferred  from  the  system  to  the  wake.  The  steady 
moment  coefficient,  Cm,  is  plotted  with  a  blue  curve.  Re=9.9T0‘*;  maximum  pitching  amplitude, 
Aa=93  deg.;  pitching  frequency,  /=5.6  Hz. 


wherein  the  force  and  torque  on  the  plate  is  significantly  augmented  by  the  formation 
and  growth  of  a  strong  LEV  on  the  suction  side  of  the  plate  (McCroskey  1982).  The  LEV 
core  continues  to  grow  in  size  and  strength  as  it  entrains  circulation  of  the  small  eddies 
generated  at  the  leading  edge.  At  the  same  time,  a  train  of  opposite-signed  vortices 
is  continuously  shed  from  the  trailing  edge — a  phenomenon  that  can  be  explained  by 
the  fact  that  the  LEV  circulation  growth  must  be  equalized  by  the  flux  of  opposite 
circulation  shed  from  the  trailing  edge  so  that  the  total  circulation  is  conserved  (i.e., 
Kelvin’s  circulation  theorem).  The  formation  of  these  discrete  eddies  has  been  commonly 
attributed  to  a  Kelvin-Helmholtz-like  instability  inherent  in  the  separated  shear  layer  by 
numerous  investigators  (e.g.  Baik  et  al.  2012;  DeVoria  &  Ringuette  2012).  The  velocity 
shear  arising  from  the  high  velocity  of  the  impinging  flow  and  the  relatively  low  velocity 
of  the  separated  free  shear  layer  at  the  sharp-edges  promotes  instability  in  the  shear  layer, 
and  therefore  leads  to  the  generation  of  the  small-scale  eddies.  Moreover,  the  shedding 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Vortex  formation  and  shedding  from  a  pitching  plate  7 

frequency  associated  with  the  small  eddies  is  much  higher  than  that  of  the  primary 
vortices  which  is  comparable  with  the  pitching  frequency. 

The  aerodynamic  torque  attains  a  maximum  at  a  «  40°  after  which  it  immediately 
begins  a  precipitous  decrease  (figure  3:  (4)).  This  torque  decrement  is  not  due  to  vortex 
leaving  the  plate,  given  that  the  attached  LEV  is  evident  at  a  ss  63°  (figure  3:  (5)). 
Rather,  we  hypothesize  that  the  saturation  of  the  LEV  circulation,  and  the  movement 
of  the  vortex  core  closer  to  the  mid-chord  point  are  responsible  for  the  decrease  in  the 
aerodynamic  torque.  This  hypothesis  will  be  revisited  in  §3.2  and  §3.4.  Despite  this  sharp 
drop  in  the  aerodynamic  torque,  the  plate  continues  to  rotate  (due  to  its  inertia)  until 
it  comes  to  rest  at  a  maximum  pitch  angle,  in  this  particular  case,  approximately  90 
degrees  to  the  free  stream  direction.  The  shear-layer  from  the  trailing  edge  begins  to  roll 
up  and  forms  a  trailing-edge  vortex  (TEV)  (figure  3:  (6)-(7)).  This  roll-up  of  the  shear 
layer  is  presumably  due  to  a  global  instability  of  the  wake  which  includes  the  influence 
of  the  LEV  passing  just  overhead.  Another  robust  feature  associated  with  the  dynamics 
of  a  strong  LEV  is  the  generation  of  small  discrete  eddies  surrounding  the  isolated  LEV 
core  in  the  wake,  as  clearly  portrayed  in  figure  3:  (7).  These  small  eddies  emerge  as 
the  vortex  filament  decays  as  a  result  of  the  Kelvin-Helmholtz  instability  (Gustafson  & 
Albert  1991).  Upon  reaching  a  maximum  angle  of  rotation,  we  observe  a  second  increase 
in  the  aerodynamic  torque  owing  to  the  development  of  a  secondary  vortex  in  the  close 
proximity  of  the  leading  edge  (figure  3:  (8)-(10)).  This  secondary  vortex  is  analogous  to 
that  observed  by  Ringuette  et  al.  (2007),  who  examined  the  LEV  formation  on  a  flat 
plate  undergoing  a  translating  motion  perpendicular  to  the  incoming  flow. 

During  the  pitch  reversal,  a  small  region  of  negative  hysteresis  is  observed  due  to  the 
fact  that  the  plate’s  return  motion  is  hampered  by  the  counter  torque  induced  by  the 
secondary  vortex.  For  applications  interested  in  energy  harvesting,  the  negative  hysteresis 
implies  that  the  energy  is  being  transferred  from  the  system  into  the  wake  and  represents 
a  drop  in  harvesting  efficiency.  We  also  note  that  the  secondary  TEV  develops  as  the 
shear  layer  emanating  from  the  trailing  edge  curls  up  in  the  wake. 

Finally,  the  elastic  restoring  toque  of  the  torsional  spring  pulls  the  plate  back  to  its 
equilibrium  angle  of  attack  (figure  3:  (11)-(13)).  Even  at  very  small  angle,  the  flow  on 
the  suction  side  of  the  plate  does  not  completely  reattach,  and  as  mentioned  eariler,  this 
results  in  a  vertical  pressure  gradient  across  the  plate  that  gives  rise  to  a  finite  torque 
at  zero  angle  of  attack.  By  symmetry,  the  pattern  repeats  during  the  second-half  of  the 
cycle. 


3.2.  Leading-edge  vortex  dynamics 

In  order  to  examine  the  effects  of  pitching  frequency  and  amplitude  on  the  ensuing  LEV 
formation,  let  us  consider  four  different  cases  of  the  reduced  frequency.  In  this  paper, 
the  reduced  frequency  is  defined  as  f*  =  7ifcAa/2Uao,  where  /  is  the  pitching  frequency 
and  Aa  is  the  maximum  pitching  amplitude  in  radians.  Note  that  cAa/2  represents  the 
leading  edge  excursion  as  the  plate  rotates  from  ao  to  amax,  during  which  the  LEV 
formation  takes  place.  This  amplitude-based  reduced  frequency  describes  the  dynamic 
similarity  between  the  leading-edge  velocity,  /cAQf/2,  and  the  characteristic  flow  velocity, 
Uoo,  which  resembles  the  inverse  of  the  advance  ratio  parameter  defined  by  Ellington 
(1984).  It  is  obvious  that  in  order  to  precisely  characterize  the  wing-tip  velocity,  one 
must  consider  both  pitching  frequency  and  stroke  amplitude — something  not  reflected  in 
the  traditional  formulation  of  the  reduced  frequency,  i.e.,  TifcjtJoo-  As  will  be  discussed 
later,  the  wing-tip  velocity  has  an  important  implication  on  the  evolution  of  the  LEV  as 
it  is  intrinsically  related  to  the  feeding  shear-layer  velocity  which  regulates  the  advection 
of  circulation  into  the  vortex  core.  The  pitching  frequency  and  amplitude  were  modulated 
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Figure  4.  (a)  Snapshots  of  the  spanwise  vorticity  fields  associated  with  various  reduced  fre¬ 
quencies  (/*=0.038  ~  0.105)  at  i?e=9.9-10^.  The  vorticity  is  scaled  by  cjUoa'-  negative,  clockwise 
(blue)  and  positive,  counter-clockwise  (red),  (b)  Corresponding  aerodynamic  torque  versus  an¬ 
gular  position  phase  portraits.  Shaded  regions  represent  the  standard  deviation  over  several 
oscillation  cycles.  All  experiments  are  performed  at  i?e=9.9T0^. 


by  tailoring  the  torsional  stiffness  and  damping  characteristics  of  the  structure.  In  each 
case,  the  time  history  of  the  angular  position  revealed  a  very  strong  sinusoidal  character 
with  a  nearly  constant  amplitude  and  a  well-defined  frequency  such  that  the  sinusoidal 
regression  curve  with  a  95%  confidence  interval  yielded  an  R-squared  value  greater  than 
0.99.  The  Reynolds  number  was  held  constant  at  9.9T0"^  in  all  experiments. 

A  strong  dependence  of  the  reduced  frequency  (defined  as  nfcjUac)  on  LEV  circulation 
and  resulting  forces  has  been  well  documented  (01  et  al.  2010;  Baik  et  al.  2012;  Granlund 
et  al.  2013;  Jantzen  et  al.  2014).  For  example,  Baik  et  al.  (2012)  and  Jantzen  et  al. 
(2014)  concluded  that  the  reduced  frequency  is  the  primary  parameter  governing  the  flow 
evolution  and  the  unsteady  force  generation  on  a  flat  plate  with  sufficient  forcing  and 
large  effective  angle  of  attack.  In  contradiction  to  their  findings,  however,  we  found  that 
the  fluid  behaviour  becomes  progressively  similar  as  f*  increases  and  nearly  independent 
of  the  reduced  frequency  for  /*>0.09,  as  can  be  inferred  from  figure  4(a).  Moreover, 
the  apparent  self-similarity  of  the  Tf  —  a  phase  trajectories  for  /*>0.07  also  reinforces 
this  assertion  (figure  4b),  which  further  implies  that  the  reduced  frequency  is  not  a 
sufficient  scaling  parameter  for  the  characterization  of  the  LEV  and  the  accompanying 
aerodynamic  torque.  However,  a  similarity  in  the  torque  trajectories  does  not  necessarily 
entail  a  similarity  in  the  force  curves,  partly  because  the  latter  is  independent  of  the 
location  of  the  centre  of  pressure.  Due  to  the  fact  that  the  pressure  gradients  acting  on 
the  plate  account  for  a  very  large  proportion  of  the  resultant  force,  as  noted  by  Albrecht 
et  al.  (2012),  the  force  should  scale  with  the  suction  pressure  induced  by  the  LEV  and 
reach  a  maximum  approximately  when  the  LEV  of  maximal  strength  is  shed  from  the 
plate,  at  which  time  the  strongest  suction  is  presumably  generated;  this  hypothesis  will 
be  verified  in  §3.4.  Because  the  LEV  circulation  increases  roughly  with  the  pitching 
frequency,  as  will  be  seen  in  §3.3,  the  force  curves  corresponding  to  /*>0.09  may  reveal 
distinct  characters  despite  the  qualitative  similarity  of  the  vorticity  fields  shown  in  figure 
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Figure  5.  (a)  Definition  sketch  of  the  vortex  centroid  with  respect  to  the  plate.  Ra'-  axial  distance 
between  the  LEV  centroid  and  the  pitch  axis;  Rn'-  shortest  distance  between  the  LEV  centroid 
and  the  axis  of  the  plate,  (b)  Evolution  of  Ra  as  a  function  of  angle  of  attack.  Different  colours 
represent  different  reduced  frequencies;  black  (/*  =0.038),  green  (/*  =0.065),  red  (/*  =0.095)  and 
blue  (/*=0.105).  (c)  Evolution  of  i?„/c  as  a  function  of  angle  of  attack.  All  experiments  are 
performed  at  i?e=9.9-10^. 


4(a).  Interestingly,  the  measurements  of  Granlund  et  al.  (2013)  have  shown  that,  for  a 
flat  plate  executing  a  linear  pitch-up  manoeuvre  from  a=0  to  90°  with  the  rotational  axis 
placed  at  3c/4,  the  force  coefficient  curves  become  similar  for  reduced  pitch  rates,  defined 
as  K=ac/2Uao,  higher  than  0.1.  As  discussed  by  the  authors,  this  observation  resonates 
with  the  unsteady  aerofoil  theory  (Leishman  2000)  which  predicts  that  the  effect  of  pitch 
rate  on  the  lifting  force  becomes  negligible  if  the  pivot  point  is  located  at  3c/4. 

To  examine  the  details  of  the  LEV  dyanmics,  the  spatial  location  of  the  LEV  centroid 
was  tracked  using  the  vortex  identification  technique  described  in  §2.  Figure  5(b)  shows 
the  evolution  of  the  normalized  axial  distance,  Rafc,  as  a  function  of  the  angle  of  attack 
(refer  to  figure  5(a)  for  notation).  For  the  lowest  frequency  case,  the  LEV  centroid  con- 
vects  more  than  c/2  aft  of  the  aerodynamic  centre,  c/4,  during  the  course  of  the  vortex 
growth  phase.  Snapshots  of  the  corresponding  spanwise  vorticity  fields  exhibit  salient 
features  of  the  separated  flow  worth  noting — streamwise  elongation  of  the  vortex  core 
and  early  separation  of  the  LEV  (figure  4(a):  top  row).  With  increasing  /*,  one  observes 
that  the  vortex  core  becomes  tighter  and  the  vortex  separation  is  delayed  to  a  higher 
angle  of  attack.  This  effect  is  clearly  illustrated  in  the  two  highest  frequency  cases  (figure 
4(a):  bottom  two  rows;  figure  5(b)),  in  which  very  tight  and  strong  LEV  cores  drift  no 
more  than  c/5  aft  of  the  aerodynamic  centre  until  the  vortices  separate  at  approximately 
aR:!60°.  This  observation  supports  the  notion  that  sufficient  rotational  acceleration  (due 
to  high  pitching  frequency)  may  contribute  to  the  stability  of  LEVs  on  wings  undergoing 
flapping  manoeuvres,  as  suggested  by  Lentink  &  Dickinson  (2009).  It  is  important  to 
note  here  that  the  qualitative  similarity  in  the  vortex  dynamics  between  these  two  cases 
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Figure  6.  Temporal  evolution  of  the  LEV  circulation  for  various  motion  profiles  at  two  different 
Reynolds  numbers.  The  circulation  and  time  are  scaled  with  (a)  the  free-stream  velocity,  Uoo 
and  (b)  the  shear-layer  feeding  velocity,  Usl- 

(despite  the  large  differences  in  /  and  Aa)  suggests  that  the  proposed  scaling  of  the 
pitching  frequency  is  appropriate  for  the  present  problem  (figure  4(a):  bottom  two  rows). 
In  addition,  the  velocity  measurements  over  a  range  of  operational  conditions  indicate 
that  the  LEV  characteristics  become  quite  similar  when  the  variation  in  f*  is  small. 

Figure  5(c)  illustrates  the  evolution  of  the  normalized  distance  between  the  LEV  cen¬ 
troid  and  the  plate  axis,  Rn/c,  as  a  function  of  the  angle  of  attack.  An  approximate 
separation  angle  can  be  inferred  from  the  figure  where  i?„/c  starts  to  increase  rapidly 
as  a  result  of  the  advection  of  the  LEV  into  the  wake.  An  important  observation,  intro¬ 
duced  earlier,  is  that  the  aerodynamic  torque  starts  to  decrease  well  before  the  vortex 
separates  from  the  plate  (figure  4a).  In  order  to  precisely  explain  this  phenomenon,  the 
relationship  between  the  LEV  formation  and  the  accompanying  aerodynamic  torque  will 
be  examined  in  §3.4  using  a  potential  flow  model. 

3.3.  Scaling  of  the  leading-edge  vortex  growth 

In  order  to  determine  the  scaling  for  the  LEV  formation  time  and  strength,  we  calcu¬ 
late  the  spanwise  circulation  of  the  LEVs  for  various  pitching  profiles.  The  experiments 
are  performed  over  a  wide  range  of  reduced  frequencies  (/*=0.038  ~  0.11),  amplitudes 
(Aa=42  deg.  ~  100  deg.)  and  at  two  different  Reynolds  numbers.  Various  LEV  circula¬ 
tion  histories  are  compared  in  figure  6(a),  where  the  circulation  and  time  are  normalized 
using  the  free-stream  velocity  and  chord.  Here,  t=0  corresponds  to  an  instant  when  the 
plate  starts  to  rise  from  the  equilibrium  angle,  a=0.  The  circulation  of  the  LEVs  increases 
monotonically  with  the  angle  of  attack  and  reaches  a  maximum  approximately  when  the 
vortex  lifts  off  from  the  plate  surface.  Subsequently  the  LEV  circulation  remains  nearly 
constant  as  the  vortex  propagates  downstream  in  the  near  wake.  This  is  a  manifestation 
of  the  fact  that  the  circulation  around  a  closed  material  contour  moving  with  the  fluid 
must  remain  constant  with  time  (i.e,  Kelvin’s  circulation  theorem).  Large  variations  in 
magnitude  and  temporal  evolution  of  the  LEV  circulation  are  clearly  evident,  in  partic¬ 
ular  between  the  low  and  high  frequency  cases,  suggesting  that  the  free-stream  velocity 
and  chord  are  not  sufficient  scaling  parameters  to  characterize  the  LEV  growth  under  a 
variety  of  conditions. 

Provided  that  the  vortex  growth  is  dictated  by  the  transport  of  circulation  flux  into 
the  vortex  core  through  the  feeding  shear-layer,  as  noted  by  Sattari  et  al.  (2012),  we 
explore  the  possibility  of  scaling  the  LEV  growth  rate  and  circulation  using  the  char¬ 
acteristic  shear-layer  feeding  velocity  at  the  leading  edge.  As  defined  in  figure  6(b),  we 
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approximate  the  shear-layer  feeding  velocity  as  the  vector  sum  of  the  local  velocity  of  the 
leading  edge  and  the  component  of  the  free-stream  in  the  direction  of  the  plate  motion: 
UsL  =  4:fcAaf2  +  Uoo  sin(a(t)).  By  construction,  this  definition  ensures  that  the  shear- 
layer  feeding  velocity  is  approximately  aligned  with  the  orientation  of  the  shear  layer  at 
the  leading  edge.  Figure  6(b)  reveals  a  good  collapse  of  circulation  histories.  The  effi¬ 
cacy  of  using  the  shear-layer  feeding  velocity  as  a  scaling  parameter  for  the  LEV  growth 
and  strength  has  been  recently  demonstrated  by  Kriegseis  et  al.  (2013)  for  a  rapidly 
accelerating  flat  plate  undergoing  plunging  motions  at  a  constant  angle  of  attack.  It 
is  worthwhile  to  mention  that  the  current  non-dimensionalization  of  vortex  circulation 
is  analogous  to  how  Dabiri  (2009)  defined  the  dimensionless  optimal  vortex  formation 
number,  T.  The  premise  of  a  universal  formation  number  is  that,  for  a  given  vortex  gener¬ 
ator,  the  maximum  vortex  strength  is  achieved  at  T«4  when  the  vortex  stops  entraining 
additional  circulation  emanating  from  the  shear  layer.  We  further  note  that  the  values 
of  vciax^iT / cU sl)  and  corresponding  dimensionless  time,  tUsh/c,  are  all  within  a  narrow 
range  of  3.5  to  4,  conforming  very  well  with  the  concept  of  universal  vortex  formation 
time  (Gharib  et  al.  1998;  Milano  &  Gharib  2005;  Dabiri  &  Gharib  2005;  Ringuette  et  al. 
2007).  Because  this  scaling  yields  an  approximately  constant  circulation  value  with  varia¬ 
tions  in  pitching  amplitude,  frequency,  and  Reynolds  number,  it  provides  a  useful  bound 
for  predicting  the  angle  of  attack,  for  a  given  pitching  frequency  and  flow  velocity,  at 
which  a  LEV  of  maximum  circulation  is  generated.  This  offers  useful  insights  as  to  how 
we  can  tailor  the  plate  kinematics  in  order  to  prolong  the  LEV  residence  time,  and  to 
leverage  the  dynamic  stall  phenomenon  to  its  fullest  advantage. 


3.4.  Potential  flow  model 

To  clarify  the  connection  between  the  vortex  evolution  and  the  resulting  aerodynamic 
torque  on  the  flat  plate,  a  potential  flow  model  with  the  addition  of  external  vortices 
(LEV  and  TEV)  is  considered  here.  The  complex  potential  function  can  be  conveniently 
derived  using  the  Joukowski  transformation,  i.e.,  z  =  (^  +  of  so  that  the  flat  plate  of 
chord  c  in  the  physical  z-plane  is  mapped  from  a  circle  of  radius  a=c/4  in  the  virtual 
(■-plane  (Batchelor  2000).  The  Gartesian  frame  is  defined  with  the  x-axis  oriented  along 
the  plate’s  chord  and  the  origin  located  at  the  mid-chord  point.  By  applying  the  circle 
theorem  (Milne-Thomson  2011),  an  expression  for  the  complex  potential  that  models 
a  stationary  flat  plate  in  the  presence  of  a  single  leading-edge  vortex  of  strength  Tlev 
located  at  (lev  and  a  trailing-edge  vortex  of  strength  Ttev  located  at  (tev  is  given  by: 


w{0 


2ni 


In 


((  (lev  \ 

C-aV((Ev/ 


-  TEV 

2m 


In 


(C  Ctev  I 

C-  aV^Ev/ 


(3.1) 


where  values  for  Tlev)  (levj  Ftev)  and  (tev  are  determined  experimentally.  In  order  to 
satisfy  the  no-penetration  condition  on  the  plate  surface,  image  vortices  with  circula¬ 
tions  — Flev  and  —Ftev  are  placed  at  the  inverse  square  points  a^/(LEv  and  a^/((Ev> 
respectively.  The  circulation  Fq  is  placed  at  the  centre  of  a  circle  in  order  to  maintain  the 
stagnation  point  at  the  trailing  edge,  whose  value  is  determined  by  imposing  the  Kutta- 
condition,  i.e.,  fizc/(i(|^=a  =  0.  This  term  merely  serves  as  a  correction  to  regularize  the 
flow  at  the  rear  edge,  and  it  should  be  distinguished  from  the  bound  circulation  of  the 
flat  plate  in  a  real  flow. 

The  aerodynamic  force  (per  unit  span)  can  be  evaluated  using  the  Blasius  formula:  — 
iFy  =  \'ipf  §c{^Ydz  (Batchelor  2000),  where  the  contour  integration  can  be  performed 
using  the  Gauchy  residue  theorem.  The  quasi-steady  approach  to  force  evaluation  can 
be  justified  based  on  the  fact  that  the  force  term  proportional  to  the  time  derivative  of 
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circulation  (i.e.,  F  (x.  dV /dt)  violates  a  number  of  physical  principles,  as  noted  by  Minotti 
(2002).  For  example,  the  author  demonstrated  that  the  temporal  change  in  the  bound 
circulation  on  a  stationary  cylinder  would  result  in  a  finite  force  that  violates  isotropy 
in  the  two-dimensional  plane.  Analogously,  in  the  case  of  a  thin  flat  plate,  a  similar 
situation  would  result  in  a  violation  of  reflection  symmetry.  Furthermore,  the  so-called 
added-mass  effect  is  expected  to  play  only  a  minor  role,  if  any,  in  the  force  generation 
owing  to  the  symmetry  of  rotational  acceleration.  For  this  reason,  the  added-mass  term 
in  the  potential  flow  model  vanishes  naturally.  The  dependence  of  the  aerodynamic  loads 
on  the  pitch  rate  or  unsteady  wing  velocity  is  indirectly  reflected  through  the  values  of 
vortex  location  and  intensity,  where,  in  §3.2  and  §3.3,  the  effects  of  wing  velocity  on  these 
quantities  are  discussed  in  detail. 

An  explicit  expression  of  the  quasi-steady  force  imparted  on  the  plate  by  the  external 
vortices  is  given  by: 


Fx  -  ^P/FleV^ 


ClBV  (  Tlbv  r^BV 

Clev  ~  \  2to(Clev  —  0?  j Clbv)  2to(Clbv  ~  Ctev) 


+ipfFi 


2ra(CLBv  -  o^/Cev)  2toCl 


St 


+ 


Clev® 


2to  (Clev-®^)^ 


Ctev  ®^  \  27Ti(Ctev  ®^/Ctev)  2re(^TEV  Clev) 


(3.2) 


-f 


2to(Ctev  ®^/Clbv)  27TzCt 


-f 


Ctev® 


2ra  (C?EV 


where  Fy  represents  the  resultant  (normal)  force  and  the  corresponding  lift  and  drag 
are  Fy  cos{a)  and  FySin(a),  respectively,  which  can  be  appropriately  normalized  by 
0.5pfU^c.  It  is  worthwhile  to  mention  that  the  last  term  in  each  of  the  curly  brack¬ 
ets  captures  what  is  known  as  Routh’s  correction  (Clements  1973),  which  rectifies  the 
velocity  discrepancy  between  the  mapping  and  physical  planes  at  the  location  of  the 
point  vortex. 

We  employed  the  Bernoulli’s  equation  to  calculate  the  pressure  distribution  around  the 
circle  in  the  C-plane,  which  is  subsequently  integrated  to  determine  the  average  location 
of  the  pressure  variation  (i.e.,  centre  of  pressure).  The  corresponding  location  on  the  flat 
plate  in  the  physical  z-plane  is  obtained  via  conformal  mapping.  The  aerodynamic  torque 
imparted  by  the  vortex  is  determined  by  taking  the  product  of  the  resultant  force  and  the 
distance  between  the  centre  of  pressure  and  the  pitch  axis.  When  the  plate  is  immersed 
in  a  real  flow,  it  is  subject  to  an  additional  torque  arising  from  the  friction  stress  and 
the  pressure  accumulation  on  the  front  surface  as  the  incoming  fluid  is  slowed  down 
and  redirected  to  turn  around  the  sharp  edges.  An  estimate  of  this  additional  torque 
is  obtained  directly  from  the  measured  steady  moment  of  the  flat  plate.  Therefore,  we 
approximate  the  net  aerodynamic  torque  as  the  sum  of  the  vortex  torque  deduced  from 
the  potential  flow  model  and  the  steady  moment  of  the  flat  plate. 

To  evaluate  the  accuracy  of  this  formulation,  we  first  consider  a  representative  case 
(corresponding  to  /*  =0.095  in  figures  4  ~  6)  to  examine  how  the  vortex  intensity  and 
stability  contribute  to  the  generation  of  the  aerodynamic  torque  on  the  plate.  Despite  a 
number  of  simplifications  and  assumptions,  the  aerodynamic  torque  evaluated  based  on 
the  inviscid  theory,  when  combined  with  the  steady  moment  of  the  flat  plate,  compares 
remarkably  well  with  that  of  the  experimental  measurement,  as  illustrated  in  hgure  7(a). 
At  large  angles  of  attack,  a  significant  portion  of  the  aerodynamic  torque  is  generated 
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{a) 


Angle  of  Attack,  a  (deg.) 


{b) 


Figure  7.  (a)  A  comparison  of  the  aerodynamic  torque  predicted  by  the  potential  flow  model 
with  the  experimental  measurement  for  the  case  f*  =  0.095  shown  in  figures  4  ~  6.  Contributions 
from  the  vortex  torque  and  the  steady  moment  are  plotted  separately,  (b)  Validation  of  the 
proposed  model  over  a  range  of  reduced  frequencies  (/*  =0.038  ~  0.11):  predictions  (solid  circles) 
and  experiments  (solid  curves). 


due  to  the  presence  of  a  strong  LEV  in  close  proximity  of  the  plate.  The  precipitous 
decrease  in  the  aerodynamic  torque  prior  to  the  vortex  separation  is  attributed  to  the 
movement  of  the  LEV  core  on  the  plate,  which  in  turn  brings  the  centre  of  pressure  closer 
to  the  rotational  axis.  The  current  analysis  suggests  that  the  contribution  of  the  TEV  on 
the  aerodynamic  loads  is  relatively  small  due  to  the  fact  that  its  circulation  is  no  more 
than  50%  of  the  maximum  LEV  circulation  and  that  the  TEV  formation  takes  place 
exclusively  in  the  wake.  However,  it  should  be  noted  that  the  inclusion  of  the  TEV  in 
the  potential  flow  model  does  seem  to  yield  a  better  estimation  of  the  centre  of  pressure, 
and  hence  improves  the  accuracy  of  the  model  at  large  angles  of  attack.  When  the  LEV 
is  very  weak  or  has  not  yet  developed  at  small  angles  (typically  less  than  15  deg.),  the 
steady  moment  constitutes  a  dominant  portion  of  the  net  aerodynamic  torque  acting  on 
the  plate. 

In  order  to  verify  the  validity  of  the  model  over  a  wide  variety  of  flow  conditions,  figure 
7(b)  compares  the  aerodynamic  torque  predicted  by  the  model  with  the  experimental 
measurements  corresponding  to  four  sets  of  data  described  in  §3.2.  It  can  be  seen  that  the 
net  aerodynamic  torque  estimated  using  our  formulation  agrees  extremely  well  with  the 
empirical  data  across  a  range  of  reduced  frequencies,  suggesting  that  the  superposition  of 
the  potential  flow  model  (which  estimates  the  torque  induced  exclusively  by  the  external 
vortices)  and  the  steady  moment  of  the  flat  plate  captures  the  essential  features  of  the 
flow  physics  associated  with  the  large-scale  oscillations.  The  significance  of  the  present 
result  lies  in  that  it  validates  the  use  of  a  quasi-steady  formulation  to  characterize  the 
unsteady  aerodynamics  of  the  flat  plate  undergoing  large-scale  oscillations.  The  two- 
dimensional  representation  of  the  vortex  exploited  in  the  model  is  justified  presumably 
because  the  flow  topology  of  the  separated  flow  structure  is  largely  two-dimensional, 
which  is  partly  ensured  in  this  study  by  the  use  of  a  high  aspect  ratio  rectangular  plate. 
We  believe  that  the  current  model  offers  an  appealing  framework  to  help  us  not  only 
elucidate  a  direct  connection  between  the  large  amplitude  oscillations  and  the  properties 
of  the  accompanying  large-scale  vortical  structures,  but  also  aid  in  the  development  of  a 
low-order  point  vortex  model  for  the  characterization  of  the  unsteady  aerodynamics  of  a 
wing. 
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4.  Concluding  remarks 

By  taking  advantage  of  the  cyber-physical  experimental  setup  and  synchronized  PIV, 
we  systematically  investigated  the  formation  time  and  strength  of  LEVs  generated  on 
a  rapidly  pitching  plate.  In  this  paper,  we  formulated  the  reduced  frequency  parameter 
based  on  the  leading  edge  excursion  during  the  pitch-up  manoeuvre.  The  velocity  mea¬ 
surements  over  a  wide  variety  of  flow  conditions  suggest  that  the  global  behaviour  of  the 
separated  flow  structures  is  strongly  controlled  by  the  reduced  frequency  for  f*  <  0.09. 
Beyond  this  threshold,  however,  the  LEV  characteristics  became  progressively  more  in¬ 
dependent  of  f* ,  from  which  we  concluded  that  the  reduced  frequency  is  not  a  sufficient 
scaling  parameter  to  characterize  the  LEV  evolution  across  a  myriad  of  flow  conditions. 

We  demonstrated  that  the  LEV  growth  and  circulation  scale  with  the  characteristic 
velocity  of  the  feeding  shear-layer  over  a  range  of  reduced  frequencies  (0.038  ~  0.11)  and 
pitching  amplitudes  (42  deg.  ~  100  deg.)  at  distinct  Reynolds  numbers.  A  connection 
between  this  scaling  and  the  concept  of  universal  vortex  formation  time  was  discussed, 
and  we  demonstrated  that  the  scaled  LEV  circulation  and  time  peak  within  a  narrow 
band  of  3.5  to  4,  revealing  an  excellent  agreement  with  the  previous  results  reported  by 
Gharib  et  al.  (1998)  and  others. 

To  elucidate  the  connection  between  the  dynamics  of  the  separated  flow  structures 
and  the  resulting  aerodynamic  torque  acting  on  the  flat  plate,  we  formulated  a  simple 
potential  flow  model  with  the  addition  of  two-dimensional  vortices.  This  model,  when 
coupled  with  the  steady  moment  of  the  flat  plate,  reproduced  the  measured  aerodynamic 
torque  remarkably  well.  Most  importantly,  we  demonstrated  that  the  precipitous  decrease 
in  the  aerodynamic  torque  occurs  before  the  point  of  LEV  separation  or  the  attainment 
of  maximum  vortex  circulation  due  to  the  movement  of  the  vortex  core  (or  centre  of 
pressure)  towards  the  rotational  axis. 

In  order  to  gain  a  deeper  understanding  of  the  LEV  dynamics,  an  assessment  of  three- 
dimensional  effects  will  be  conducted  in  the  future.  On  a  different  level,  we  will  attempt 
to  optimize  the  plate  kinematics  with  the  aim  of  prolonging  the  LEV  residence  time  and 
augmenting  the  vortex  strength  by  implementing  appropriate  structural  nonlinearities  in 
the  system  and  exploring  the  effects  of  pivot  location. 
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